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A  phase  field  theory  for  coupled  twinning  and  fracture  in  single  crystal  domains 
is  developed.  Distinct  order  parameters  denote  twinned  and  fractured  domains, 
finite  strains  are  addressed  and  elastic  nonlinearity  is  included  via  a  neo-Hookean 
strain  energy  potential.  The  governing  equations  and  boundary  conditions  are 
derived;  an  incremental  energy  minimization  approach  is  advocated  for  prediction 
of  equilibrium  microstructural  morphologies  under  quasi-static  loading  protocols. 
Aspects  of  the  theory  are  analysed  in  detail  for  a  material  element  undergoing 
simple  shear  deformation.  Exact  analytical  and/or  one-dimensional  numerical 
solutions  are  obtained  in  dimensionless  form  for  stress  states,  stability  criteria  and 
order  parameter  profiles  at  localized  fractures  or  twinning  zones.  For  sufficient 
applied  strain,  the  relative  likelihood  of  localized  twinning  vs.  localized  fracture 
is  found  to  depend  only  on  the  ratio  of  twin  boundary  surface  energy  to  fracture 
surface  energy.  Predicted  criteria  for  shear  stress-driven  fracture  or  twinning  are 
often  found  to  be  in  closer  agreement  with  test  data  for  several  types  of  real 
crystals  than  those  based  on  the  concept  of  theoretical  strength. 
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1.  Introduction 

Cleavage  fracture  and  deformation  twinning  are  two  fundamental  inelastic  deformation 
mechanisms  that  may  occur  in  crystals  when  stressed  above  their  elastic  limit.  Cleavage 
fracture  is  defined  here  as  fracture  on  specific  crystallographic  planes,  typically  those  of 
lowest  surface  energy  [1].  Deformation  twinning,  also  known  as  mechanical  twinning,  is 
defined  here  as  twinning  induced  by  mechanical  stress  [2,3].  Both  of  these  anisotropic 
mechanisms  involve  deformation  on  specific  planes  (the  cleavage  plane  for  fracture  or  the 
habit  plane  for  twinning)  and  tend  to  be  dissipative  or  thermodynamically  irreversible, 
though  reversible  or  “elastic”  twinning  is  possible  in  some  materials  [4].  Local  stress 
concentrations  due  to  fracture  (e.g.  at  crack  tips)  can  induce  twinning,  and  vice  versa 
[2,5-7].  Constitutive  models  accounting  for  both  phenomena  are  needed  for  the  prediction 
of  the  onset  of  inelasticity  and  subsequent  mechanical  response  in  materials  of  interest  for 
engineering  applications,  including  many  metals  and  ceramics. 
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Numerous  analytical  and  computational  models  have  been  posited  to  describe  fracture 
or  twinning,  with  resolutions  ranging  from  discrete  atoms  to  macroscopic  continua.  In  the 
context  of  continuum  mechanical  theories,  phase  field  models  demonstrate  a  number  of  use¬ 
ful  characteristics:  structural  transformations  occur  naturally  via  energy  minimization  [8]  or 
Ginzburg-Landau  type  kinetics  [9],  without  the  need  for  ad  hoc  evolution  criteria;  relatively 
few  material  parameters/properties  are  required;  and  numerical  solutions  are  regularized 
(i.e.  are  rendered  mesh-size  independent)  by  a  user-defined  length  parameter  associated 
with  interfacial  width.  Certain  regularized  variational  models  incorporating  linear  elasticity 
may  also  demonstrate  so-called  gamma  convergence  towards  the  Griffith  theory  of  fracture 
as  the  regularization  zone  shrinks  to  a  singular  surface  [10,11];  however,  such  convergence 
has  apparently  not  been  completely  proven  for  nonlinear  elasticity  [12].  As  detailed  in  the 
following  paragraphs,  phase  field  theories  have  already  been  developed  and  implemented  for 
independent  modelling  and  simulation  of  fracture  and  deformation  twinning,  but  no  model 
seems  to  have  been  reported  previously  that  simultaneously  addresses  both  phenomena,  in 
the  same  analysis,  with  distinct  order  parameters  for  twinned  and  fractured  domains.  The 
development  of  such  a  coupled  phase  field  theory  is  the  objective  of  the  present  work. 

Phase  field  theories  for  fracture  [13-24],  or  variational  gradient-type  theories  based  on 
energy  minimization  similar  to  phase  field  approaches  [10-12,25],  have  been  under  contin¬ 
uous  development  since  the  early  2000s.  Models  have  been  developed  with  distinct  order 
parameters  for  fracture  and  electric  polarization  [26-29],  and  the  latter  can  be  physically 
related  to  formation  of  twinned  domains  in  ferroelectrics;  this  sort  of  electromechanical 
twinning  differs  from  deformation  twinning  associated  with  partial  dislocations  considered 
herein.  Though  most  of  the  aforementioned  fracture  theories  are  posited  in  the  context  of 
small  strains,  several  address  finite  deformations  and  nonlinear  elasticity  [12,15],  aspects 
which  are  also  addressed  by  the  theory  developed  in  the  present  paper.  Also  newly  addressed 
in  this  work  is  the  issue  of  stability  of  the  nonlinear  phase  field  fracture  theory,  extending 
linear  elastic  analysis  of  tension  considered  in  prior  work  [14,30-35]  to  nonlinear  analysis 
of  finite  simple  shear. 

Though  not  as  large  as  the  literature  on  phase  field  fracture  models,  phase  field  models 
of  deformation  twinning  have  also  been  documented,  with  some  theories  based  on  geo¬ 
metrically  linear  elastic  response  [36-38].  Geometrically  nonlinear  phase  field  models  of 
twinning  include  those  of  the  present  authors  [5,8,39]  applied  to  a  single  twin  system  and 
[40]  applied  to  martensitic  transformations. 

Developed  in  Section  2  of  the  present  paper  seems  to  be  the  first  phase  field  theory 
accounting  for  both  fracture  and  deformation  twinning  wherein  each  mechanism  is  repre¬ 
sented  by  a  distinct-order  parameter.  The  model  is  variational  or  quasi-static  in  nature,  and  is 
intended  for  eventual  use  in  finite  element  simulations  that  seek  equilibrium  solutions  using 
incremental  energy  minimization  of  a  total  energy  functional  that  accounts  for  elastic  strain 
energy  and  surface  energies  of  cleaved  material  and  twin  boundaries.  The  general  theory 
combines  aspects  of  previous  individual  nonlinear  models  for  twinning  [8,39]  and  fracture 
[15].  The  minimal  necessary  material  parameters  are  as  follows:  the  twinning  habit  plane 
and  twinning  eigen-shear,  the  twin  boundary  surface  energy,  the  crack  plane  normal  (for 
anisotropic  fracture),  the  surface  energy  of  fracture  and  the  elastic  constants.  A  numerical 
parameter  is  also  needed  for  regularization  such  that  the  minimum  twin  boundary/crack 
width  is  finite.  Via  appropriate  choice  of  energy  functional,  the  model  accounts  for  possible 
coupling  between  fracture  and  twin  evolution. 
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The  remaining  sections  of  the  paper  consider  application  of  various  aspects  of  the  general 
theory  of  Section  2  towards  a  one-dimensional  simple  shear  problem.  The  analysis  assumes, 
a  priori,  that  twinning  and/or  fracture  are  confined  to  occur  in  the  direction  of  applied 
shear  deformation  on  a  single  plane  upon  which  the  resolved  shear  stress  is  maximum 
(prior  to  nucleation)  or  nearly  maximum  (after  nucleation).  Only  a  single  twin  system  - 
with  its  habit  plane  aligned  with  the  plane  of  maximum  shear  stress  -  and  a  single  mode 
II  crack  are  considered.  From  a  pragmatic  mathematical  perspective,  these  assumptions 
permit  a  detailed  and  thorough  analysis  of  the  nonlinear  elasticity  problem  that  would  not  be 
possible  for  more  complex  two-  and  three-dimensional  deformation  states  involving  activity 
of  many  twinning  systems  and/or  many  cleavage  planes.  Real  crystalline  materials  tend  to 
have  multiple  available  twinning  systems  or  twin  variants  [2],  but  the  present  restriction 
to  a  single  system  is  physically  justified  for  the  simple  shear  problem  since  deformation 
twinning  is  typically  driven  by  resolved  shear  stress  [3,8,41].  For  these  boundary  conditions, 
the  resolved  shear  stresses  on  other  twin  systems  in  the  material  would  not  be  sufficiently 
large  for  simultaneous  nucleation  of  multiple  systems;  similarly,  twinning  in  calcite  loaded 
in  confined  direct  shear  [42]  or  indentation  [4]  demonstrates  activity  of  only  a  single  system 
associated  with  maximum  shear  stress  despite  the  multiplicity  of  systems  available  in  the 
crystal.  For  the  above  reasons  and  for  conciseness  of  notation,  the  theory  of  Section  2  only 
accounts  for  one  potential  twin  system;  to  address  more  general  problems,  extension  of  the 
theory  to  account  for  multiple  twin  systems  described  by  distinct-order  parameters  can  be 
achieved  following  methods  in  [8,36-38]. 

Regarding  fracture,  most  brittle  materials  tend  to  fail  initially  in  practice  by  extensional 
(mode  I)  cracking  [43],  with  crack  faces  oriented  for  opening  along  directions  of  maximum 
tensile  principal  stress.  For  a  state  of  pure  shear  stress,  the  maximum  tensile  stress  direction 
is  oriented  at  an  angle  of  45°  from  the  shearing  direction,  i.e.  from  the  direction  presumed 
a  priori  in  the  mode  II  analysis  considered  later  in  this  work,  with  the  magnitude  of  the 
maximum  tensile  stress  equal  (or  nearly  so)  to  the  applied/resolved  shear  stress.  However, 
the  shearing  mode  of  fracture  would  be  preferred  in  a  material  with  a  lower  fracture  strength 
for  tangential/sliding  modes  than  for  tensile  modes,  or  for  a  highly  anisotropic  (e.g.  mica¬ 
like)  solid  oriented  with  cleavage  planes  parallel  to  the  shearing  direction.  Furthermore, 
the  study  of  pure  shearing-type  fracture  is  of  theoretical  interest  because  results  can  be 
compared  with  other  models  of  “theoretical  shear  strength”  as  originated  by  Frenkel  [44] 
and  used  to  provide  fundamental  insight  into  aspects  of  atomic  bonding  [45]  and  upper 
bounds  for  stresses  leading  to  homogeneous  dislocation  nucleation  and  glide  [46].  Values 
of  theoretical  shear  strength  obtained  under  conditions  of  pure  or  simple  shear  loading  can 
be  used  in  failure  criteria  for  brittle  materials  undergoing  other  loading  protocols  in  which 
all  three  principal  stresses  are  compressive,  wherein  a  maximum  tensile  stress  criterion 
would  be  ineffective.  For  example,  theoretical  shear  strength  arguments  have  been  used 
to  describe  the  onset  of  fracture  in  shock  compression  [47,48]  (and  see  later  Section  4  of 
this  work),  and  shear  fracture  strength  is  of  high  importance  for  indentation  and  ballistic 
loading  [49]  in  regions  of  material  far  from  free  surfaces  where  tensile  cracking  may  occur, 
as  well  as  in  special  direct  shear  experiments  [42]  where  confinement  precludes  opening 
of  mode  I  cracks.  The  analytical  solutions  derived  later  do  not  explicitly  account  for  large 
compressive  hydrostatic  pressures  that  accompany  planar  shock  compression,  for  example, 
but  theoretical  shear  strengths  derived  from  a  tractable  analysis  of  the  simple  shear  problem 
may  be  transferable,  as  has  been  assumed  elsewhere  [47,50],  to  the  former  type  of  (shock) 
loading  that  cannot  be  subjected  to  a  complete  analytical  treatment  using  the  phase  field 
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approach.  Finally,  analytical  solutions  to  the  simple  shear  problem  with  a  single  pure  mode 
II  crack  enable  isolation  of  deviatoric  features  of  the  model  that  can  be  used  for  validation 
of  numerical  methods  (e.g.  finite  element  solutions  to  be  obtained  in  future  work)  and 
comparison  with  sharp  interface  (e.g.  Griffith-type,  see  Section  4  representations;  similar 
studies  have  been  undertaken  for  shear  in  a  geometrically  linear  setting  [25]. 

Interactions  between  fractured  and  transformed/twinned  zones  have  been  reported  for 
crystalline  solids  undergoing  simultaneous  fracture  and  twinning  or  phase  transformations, 
and  such  interactions  have  been  modelled  elsewhere  via  phase  field  numerical  simulations  of 
ferroelectrics  [29]  and  austenitic-martensitic  steels  [51].  Complex  stress  states  may  emerge 
as  a  result,  and  transformation  kinetics  influence  crack  paths  and  vice  versa.  The  present 
complete  theory,  if  implemented  numerically  for  two-  or  three-dimensional  simulations, 
would  allow  for  prediction  of  incremental  equilibrium  states  corresponding  to  such  phe¬ 
nomena,  including  crack  branching  affected  by  twinning  under  quasi-static  loading,  for 
example.  However,  the  particular  analytical  solutions  derived  for  simple  shear  with  cracks 
and  twins  aligned  parallel  to  a  single  direction  do  not  permit  consideration  of  such  complex 
stress  states  and  changes  in  crack  paths. 

The  governing  equations  and  nondimensional  forms  are  derived  for  the  simple  shear 
problem  admitting  both  fracture  and  twinning  in  Section  3.1.  The  possibility  of  damage 
and  twinning  to  appear  sequentially  or  simultaneously  and  the  relative  likelihood  of  com¬ 
plete  mode  II  fracture  vs.  localized  twinning  are  addressed  in  Section  3.2.  Section  3.3 
considers  mode  II  fracture  in  the  absence  of  twinning:  nondimensional  analytical  solutions 
for  homogeneous  damage,  stability  criteria  and  localized  complete  fracture,  as  well  as 
numerical  solutions  for  inhomogeneous  damage.  Twinning  under  simple  shear  in  the  absence 
of  damage  is  then  considered  in  Section  3.4  as  another  particular  case:  nondimensional 
analytical  solutions  for  homogeneous  elasticity,  stability  and  numerical  solutions  for  local¬ 
ized  complete  twinning.  Predicted  critical  stresses  associated  with  shear  fracture  or  twin 
nucleation  are  compared  with  other  models  and  data  from  the  literature  for  several  real 
crystalline  materials  in  Section  4.  Conclusions  follow  in  Section  5. 

Notation  of  modem  continuum  mechanics/physics  is  used.  Vectors  and  tensors  are 
written  in  bold  font  and  are  referred  to  a  Cartesian  coordinate  system.  Superscripts  T, 
—  1  and  — T  denote  the  transpose,  inverse  and  inverse-transpose.  In  indicial  notation  with 
summation  over  repeated  indices,  the  dot  product,  double-dot  product  and  outer  product 
obey  a  -  b  —  axbx,  A  :  B  —  AjjBjj,  and  (a  ®b)[j  —  ajbj. 


2.  Nonlinear  phase  field  theory 

The  present  theory  essentially  combines  previous  phase  field  approaches  for  twinning  [8,39] 
and  fracture  [15]  in  a  way  general  enough  to  permit  physically  realistic  coupling  of  the  two 
phenomena.  A  geometrically  nonlinear  theory  is  considered,  accounting  for  finite  strains, 
finite  rotations  and  nonlinear  elastic  response. 


2.1.  Order  parameters  and  kinematics 

Let  X  denote  material  coordinates  and  t  time.  Let  r}{X,t)  denote  the  order  parameter 
associated  with  twinning  on  a  single  twin  system  (see  e.g.  [8]  for  extension  to  multiple  twin 
systems): 
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ri(X,t)  —  OVZ  e  parent, 
ri{X,  t)  —  IVZ  e  twin, 

ri{X,  t)  e  (0,  1)VZ  e  twin  boundary.  (2.1) 

Let  f(Z,f)  denote  the  order  parameter  associated  with  fracture  (possibly,  but  not 
necessarily,  concentrated  on  a  single  cleavage  plane): 

f  (Z,  f)  =  OVZ  e  intact  solid, 

§(Z,  t)  =  IVZ  e  failed  domain, 

§(Z,  t)  e  (0,  1)VZ  e  cohesive  boundary.  (2.2) 

Motion  is  Jc(Z,  f)  =  Z  +  h(Z,  t)  with  u  the  displacement.  The  deformation  gradient 
is  the  two-point  tensor 

F  =  Vx  =  1  +  Vm,  (2.3) 

where  V/f(-)  =  3(-)/9Z/f  is  the  material  gradient.  The  deformation  gradient  is  split 
multiplicatively  as  [8] 

F  =  (2.4) 

Twinning  deformation  is  the  simple  shear 

F’^(r})  =  1  4-  [(p{r\)YQ\s  (g)  m,  det  F''  =  1  4-  cpyos  ■  m  —  1.  (2.5) 

Orthogonal  unit  vectors  (in  material  coordinates)  in  the  twinning  direction  and  normal 
to  twinning  plane  are  s  and  m,  the  magnitude  of  maximum  twinning  shear  is  yo,  and 
q)  :  [0,  1]  — ^  [0,  1]  is  a  continuous  and  monotonically  increasing  interpolation  function 
satisfying  ^(0)  =  0  and  ^(1)  =  1.  Define  the  symmetric  elastic  material  deformation 
tensor  and  the  Jacobian  determinant: 

C  =  F^^F^,  J  =  det  F^  =  det  F  =  VdetC.  (2.6) 

An  incremental  or  quasi-static  theory  is  the  focus,  with  explicit  dependence  of  field 
variables  on  t  omitted. 

A  few  remarks  regarding  compatibility  are  in  order.  Neither  F^  nor  F'*  is  generally 
integrable  to  a  vector  field;  i.e.  elastic  and  twin  deformation  mappings  in  that  case  are 
generally  anholonomic  [52,53],  meaning 

VxFVO,  VxF'^-^^O,  (2.7) 

where  V  is  the  gradient  with  respect  to  spatial  coordinates  [i.e.  V^.(•)  =  d{-)/dxk  — 
V7(  )FJ^'].  Total  deformation  gradient  F  and  its  inverse  are  compatible  by  definition: 

VxF  =  VxVx  =  0,  V  X  F“' =  V  X  VZ  =  0.  (2.8) 

In  contrast  to  fracture  mechanics  models  wherein  discrete  cracks  are  resolved  explicitly 
in  terms  of  jump  discontinuities  in  displacement  (leading  to  singularities  in  F),  in  the  phase 
field  approach  the  displacement  field  h(Z)  is  continuous  and  continuously  differentiable. 
As  explained  later  in  Section  2.3,  failure  is  captured  via  a  reduction  in  local  elastic  stiffness 
(i.e.  an  increase  in  compliance)  at  a  given  material  point  Z,  similar  to  continuum  damage 
mechanics  theories.  Accordingly,  an  advantage  of  the  phase  field  method  over  computational 
models  involving  discontinuous  fields  (e.g.  cohesive  finite  elements)  is  that  conventional 
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finite  element  or  finite  difference  schemes  can  be  invoked  in  the  former,  requiring  no 
special  treatment  or  tracking  of  surfaces  of  discontinuity.  The  present  phase  field  theory  also 
differs  from  other  finite  strain  continuum  damage  mechanics  models  wherein  incompatible 
components  of  the  deformation  gradient  are  used  to  represent  irreversible  deformation 
associated  with  local  fractures  within  a  volume  element  of  material  [54—56]. 


2.2.  Governing  equations 

The  total  energy  functional  for  an  initially  homogeneous  body  with  uniform  material 
properties  is 

P,  t)  =  f  [W(Vx,  q,  I)  +  fir,,  Vr,,  Vt)]dy.  (2.9) 

Jq 

In  the  reference  configuration,  the  body  occupies  volume  V  and  is  enclosed  by  surface 
S  with  unit  outward  normal  n.  In  the  nonlinear  theory,  strain  energy  per  unit  reference 
configuration  volume  W  is  of  the  form 

W  ^W[CiF,r,),r,,^].  (2.10) 


Stationary  points  of  the  total  energy  functional  4'  are  obtained  as  solutions  of 


where,  from  (2.9), 


S^r  —  (t  (t  ■  Sx  +  rSr,  +  sS^)dS, 

hn 


(2.11) 


54/  = 


l\ 

Jn  L 


dW 

Jf 


dW  dW  df  df  df  9/ 

:  VSx  + - Sr,  + - +  —Sr,  +  —S^  +  ■  WSr,  + 

dr,  dr,  '  S$  ^  dVr,  9Vf 


dV. 

(2.12) 


Assuming  fields  are  sufficiently  smooth,  the  following  general  local  equilibrium  equations 
(i.e.  Euler-Lagrange  equations)  and  boundary  conditions  are  then  derived  by  extending,  for 
example,  the  procedure  involving  the  divergence  theorem  and  integration  by  parts  detailed 
in  [8]: 


V  ■ 


dW 

Jf 


dW 

=  VP=0,  P  =  - 

>74  dF 


>74 


9/  9/ 

—  -  V  ■  — ^  + 
dr,  dVr, 

t  —  P  ■  n,  r  — 


dW 

dr, 

9/ 

dVr, 


df  df  dW 

=  0,  ^  -  V  ■  H - 

av^  at 

9/ 

■  «,  5  =  - ■  n. 

av^ 


F,r, 


=  0; 


(2.13) 

(2.14) 

(2.15) 


The  traction  per  unit  reference  area  is  t\  the  conjugate  force  to  order  parameter  r,  is  r,  the 
conjugate  force  to  order  parameter  f  is  s.  These  equations  apply  for  the  general  framework 
of  (2.9),  with  (2.13)  and  (2.14)  holding  within  domain  and  (2.15)  holding  on  its  external 
boundary  3^2. 

The  interfacial  energy  per  unit  reference  volume  considered  here  is  of  the  general  form 


fir,,  f ,  Vp,  V§)  =  Mr,,  ^)  +  go(t)  +  fiid,  Vt?)  +  ^1(77,  Vf).  (2.16) 
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Taking  this  into  account,  equilibrium  Equation  (2.14)  and  boundary  conditions  (2.15)  can 
be  expressed  as 


dr]  dr] 
9/o  9go 
9?  9? 

t  —  P  ■  n, 


-  V  ■ 


9gi 

dr] 

9/i  9gi 
9/i 


9/i 

dVi] 


-  V 


dVi] 


■  n, 


dW 

dr] 

9gi 

9V^ 

9gi 


=  0; 
F.^ 

dW 


9^ 


FJI 


9Vt 


•  n. 


=  0; 


(2.17) 

(2.18) 
(2.19) 


2.3.  Compressible  neo-Hookean  elastieity 

Similar  to  previous  phase  field  models  for  twinning  [5,39],  compressible  neo-Hookean 
elasticity  is  now  considered,  where  for  modelling  elasticity  and  fracture, 

1T=  i[A.(ln/)2-/r(21n7-trC  +  3)].  (2.20) 

Lame  coefficients  /r  and  X  depend  on  §  and  possibly  J: 

/r(|)  =  Mo{C  +  (l-O[l-0(f)]},  k  =  k(^, /)  -  |/x(§).  (2.21) 

Here,  cp  :  [0,  1]  — [0,  1]  is  a  continuous  monotonic  interpolation  function  satisfying 
(piO)  —  0  and  (/>(1)  =  1,  and  (;■  is  a  small  scalar  constant  subject  to  0  <  C  «  1  .  A  particular 
example  function  is 

</.(§)  =  1 -(1-1)2.  (2.22) 

The  shear  modulus  p.  degrades  from  its  reference  value  /tq  to  a  minimum  value  C,pQ  in 
fully  fractured  domains.  The  constant  C,  ensures  that  the  material  maintains  a  minimal 
residual  stiffness  when  fully  fractured  [22,26];  for  solution  of  some  (but  not  all)  boundary 
value  problems,  it  may  be  needed  for  numerical  stability  and  is  not  considered  a  material 
property.  For  conciseness,  a  degradation  function  0  is  introduced: 

0(§)  =  C  +  (1  -  0[1  -  0(f)]  ^  M  =  MO0.  (2.23) 

Note  that  the  prescription  of  a  regularized  variational  formulation  of  fracture  in  [10]  is 
recovered  when  the  residual  stiffness  is  omitted: 

c  =  0  ^  0(?)  =  (1  -  ?)2  ^  =  [1  -  0(t)]/ro  =  (1  -  f )Vo.  (2.24) 

A  quadratic  degradation  of  elastic  stiffness  with  damage  has  likewise  been  used  in  a  number 
of  other  phase  field  and  variational  models  of  fracture  [11,14,20,22].  The  bulk  modulus  k 
degrades  in  tension  but  not  in  compression,  in  order  to  prohibit  interpenetration  [25]: 

k  =  (ko  +  I /ro)[{7  -  l)0(f )  +  (1  -  7)*],  (2.25) 

where  ko  is  the  constant  initial  Lame  modulus  and  (x)  =  IVx  >  0,  (x)  =  OVx  <  0, 
(x)*  =  IVx  >  0,  and  (x)*  =  OVx  <  0.  The  above  neo-Hookean  model  is  restricted 
to  isotropic  elasticity,  a  notable  approximation  for  real  crystals.  For  anisotropic  elastic 
constants,  the  model  must  be  extended  to  account  for  transformation  (reflection  or  rotation) 
of  the  reference  frame  of  the  crystal  lattice  commensurate  with  twinning,  as  has  been 
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elaborated  for  a  nonlinear  model  based  on  the  Green  strain  [8] .  However,  since  Green  strain- 
type  theories  incorporating  higher  order  elastic  constants  tend  towards  instability  under  large 
compression  [57,58],  formulation  of  a  nonlinear  anisotropic  theory  incorporating  another 
strain  tensor  such  as  the  material  Eulerian  strain  [59,60]  or  material  logarithmic  strain  [61] 
may  be  prudent  in  the  future. 

Further  justification  of  the  use  of  a  nonlinear  elastic  model,  and  compressible  neo- 
Hookean  elasticity  in  particular,  is  now  discussed.  In  brittle  solids  undergoing  fracture  (but 
not  twinning),  linear  elasticity  may  be  sufficient  if  failure  occurs  before  large  deformations 
are  attained;  recent  phase  field  simulations  [15]  have  demonstrated  very  similar  results  for 
crack  nucleation  in  linear  and  nonlinear  elastic  brittle  solids.  However,  the  present  work 
is  concerned  with  both  fracture  and  deformation  twinning,  and  accurate  depiction  of  the 

latter  usually  requires  a  finite  deformation  description  [62].  Twinning  shear  deformation  yo 

fi 

is  often  large,  for  example  a  value  of  for  the  most  common  mode  in  cubic  crystals  [2]. 
Previous  phase  field  simulations  [5,39]  have  demonstrated  differences  in  predictions  for  twin 
nucleation  and  growth  among  linear  elastic  and  neo-Hookean  models.  The  neo-Hookean 
model  used  herein  has  also  been  shown  to  provide  a  reasonable  two-parameter  depiction 
of  the  increasing  bulk  stiffness  of  some  crystalline  solids  (calcite  and  sapphire)  under  finite 
volumetric  compression,  in  contrast  to  linear  elasticity  wherein  the  bulk  modulus  is  constant. 
Neo-Hookean  elasticity  has  also  been  incorporated  elsewhere  in  a  nonlinear  phase  field 
model  of  twinning  in  crystalline  CuAlNi  (martensitic  phase)  [40].  The  distinction  between 
the  compressible  neo-Hookean  potential  used  here  and  incompressible  neo-Hookean  elas¬ 
ticity  often  used  for  polymers  such  as  plastics  and  rubber-type  materials  is  also  noted.  In 
ductile  crystalline  solids,  crack  tip  plasticity  often  ensues  prior  to  attainment  of  large  elastic 
deformations,  rendering  consideration  of  dislocation  motion  more  important  than  nonlinear 
elastic  effects.  As  remarked  later  in  Section  4,  extensions  of  the  current  model  framework 
would  be  necessary  for  modelling  dislocation  plasticity  distinct  from  twinning. 

For  twinning  on  a  single  plane,  with  a  double-well  potential  and  isotropic  surface  energy, 

(2.26) 

/i(t,V/7)=/r(|):  (VpOVp),  *:(§)  =  (2.27) 

where  the  constants  A  and  kq  are  related  to  the  twin  boundary  surface  energy  F  and  the 
thickness  I  as  [8] 

A  =  12r//,  Ko  =  |r/.  (2.28) 

Note  that  /o  and  k  can  degrade  with  f  via  function  t;  this  can  prevent  a  twin  boundary  that 
crosses  a  free  surface  from  having  duplicate  surface  energy,  for  example.  The  degradation 
function  t :  [0,  1]  — [0,  1]  obeys 

((0)  =  1,  r  =  d?/dt.  (2.29) 

A  candidate  choice  is  t  =  0,  meaning  twin  boundary  energy  and  elastic  shear  modulus 
degrade  analogously. 

For  cleavage  fracture  on  a  preferred  plane,  it  is  assumed  that  the  orientation  of  such  a 
plane  is  known  a  priori.  For  example,  this  would  be  a  plane  or  family  planes  of  low  surface 
energy  or  low  intrinsic  strength  in  a  crystal.  Considered  herein  is  a  single  plane.  A  unit 
vector  (in  material  coordinates)  normal  to  this  cleavage  plane  is  M.  Let  B,  coq  and  denote 
material  constants.  For  the  present  fracture  model. 
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goi^)  =  Bf,  (2.30) 

gi(V§)  =0):  (V^(g)  V^),  <y  =  «o[l  +  /0(l-M0M)].  (2.31) 


The  crack  thickness  h  and  fracture  surface  energy  T  are  related  as  [22] 

B  =  0)0  =  Th.  (2.32) 


Setting  P  ^  \  greatly  increases  the  surface  energy  of  fracture  on  planes  not  normal  to 
M;  /S  is  thus  a  penalty  factor,  not  a  typical  material  property.  Setting  /j  =  0  results  in 
isotropic  damage.  This  isotropic  continuum  representation,  with  go  oc  has  been  shown 
to  converge  to  the  correct  surface  energy  of  a  singular  surface  as  1  ^  0  [10,11]  and  has 
been  used  in  a  number  of  recent  continuum  models  [12,14,20,22,25]. 

Consider  now  the  specific  constitutive  model  of  (2.20)-(2.32)  incorporating  com¬ 
pressible  neo-Hookean  elasticity,  a  double-well  potential  for  twinning  with  eigen-strain 
interpolator  q){rj)  and  a  quadratic  potential  for  fracture  with  degradation  function  0(^). 
Euler-Lagrange  equations  for  order  parameters  (2.17)  and  (2.18)  become,  with  t  and  a 
interpreted  as  elastic  driving  forces  associated  with  twinning  and  fracture. 


r  = 

?  = 


dW 

dr] 

dW 

W 


=  2ko{?($)V^>]  +  ■  Vr]}  -  2Arj(l  -  3/?  +  2,]^)c($); 

=  2aJo[V^^  +  /3(V^$  -  M  0  M  :  VV  $)] 
-2B$-[>colVr]l^  +  Ar]^(l-r]f]r($). 


(2.33) 


(2.34) 


Primes  denote  derivatives  of  functions  of  one  variable,  e.g.  (p'  =  d(p/drj  and  (p"  =  d^(p/dr]^. 
For  simple  problems  (e.g.  one-dimensional  problems  such  as  those  considered  later  in 
Section  3),  it  may  be  possible  to  find  solutions  of  (2.13)  and  (2.33)-(2.34)  directly  and 
analytically.  General  two-  and  three-dimensional  boundary  value  problems  will  require 
numerical  solutions,  which  can  be  found  incrementally  using  the  finite  element  method  [8] 
wherein  boundary  conditions  for  displacement  and  order  parameters  (or  their  conjugate 
forces)  are  updated  during  each  load  increment.  For  certain  boundary  conditions,  candidate 
solutions  [i.e.  fields  x(X),  r]{X),  ^(X)]  are  obtained  that  minimize  'k(x,  rj,  f)  and  thus 
satisfy  (2.11).  However,  because  4'  is  generally  nonconvex,  multiple  (local)  minima  may 
exist,  and  therefore  solutions  may  be  nonunique;  when  4'  of  such  a  local  minima  exceeds 
the  global  minimum  energy,  then  such  a  solution  is  said  to  be  metastable.  To  ensure 
irreversibility  of  crack  extension,  constraints  such  as  S^(X)  >  0  for  |(Z)  exceeding  some 
threshold  value  can  be  imposed  [12,15]. 

A  number  of  different  interpolation  functions  <p  have  been  used  for  interpolation  of 
twinning  shear  [5,8,9,40].  A  particular  function  considered  later  is 

cp(,])  =  r]H3  -  2r]),  (2.35) 


with  derivatives 

(p'{r])  ^  6r](l  -  r]),  (p" {-q)  ^  6{l  -  2r]) ,  (2.36) 

and  useful  properties 

<p{0)  =  0,  ¥)(!)  =  1,  = 


^'(0)  =  <p\\)  =  0;  cp"{\)  =  0.  (2.37) 
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For  completeness,  derivatives  and  special  values  of  function  0(^)  of  (2.23)  with  0  of  (2.22) 
are  listed: 

0'(^)  =  -2(l-O(l-?),  0"(§)  =  2(1  -  C);  (2.38) 

0(0)  =  1,  0(1)  =  ^;  0'(O)  = -2(1  -  C),  0'(1)  =  O.  (2.39) 

Further  remarks  on  a  number  of  scalar  functions  entering  the  total  energy  functional  are 
now  given,  specifically  /o,  f\,  go,  gi,  0,  0,  t  and  (p.  None  of  these  functions  is  determined 
uniquely  by  mathematical  or  physical  criteria,  though  relatively  simple  forms  for  all  are 
chosen  that  satisfy  stated  necessary  conditions  associated  with  endpoint  values,  continuity 
and  smoothness,  and  possible  symmetry.  In  particular,  /o  includes  a  double-well  potential 
typically  used  in  phase  field  theory  [8,39] ,  though  other  higher  order  polynomial  forms  have 
been  prescribed  elsewhere  for  twinning  energetics  [36,38].  Quadratic  forms  prescribed  for 
gradient  terms  f\  and  are  standard  in  phase  field  models.  Quadratic  term  go  is  used 
in  conjunction  with  degradation  function  0  because,  for  linear  elasticity  in  the  absence 
of  twinning,  these  provide  gamma  convergence  of  the  regularized  fracture  energy  to  the 
surface  energy  in  discrete  (Griffith-type)  fracture  mechanics  [10,11].  The  model  used  here 
for  stiffness  degradation  in  tension  and  shear,  incorporating  0  with  residual  strength  ^ ,  agrees 
with  representations  used  elsewhere  [12,15];  however,  it  has  been  verified  in  variational 
simulations  of  fracture  [25]  that  quantitative  predictions  depend  on  choice  of  Q  and  thus 
implicitly  on  choice  of  0.  Coupling  function  t  is  new  to  the  present  theory  (since  no  other 
coupled  fracture-twinning  phase  field  theory  is  known  to  the  authors)  and  its  choice  does 
influence  the  availability  of  analytical  solutions,  as  discussed  in  Section  3.2.  Interpolation 
function  (p  has  been  used  in  other  finite  strain  models  involving  twinning  [5,8,39]  and  phase 
transformations  [9];  prior  work  has  verified  that  use  of  a  different  interpolation  function 
(specifically  an  exponential  function  of  Fermi-Dirac  type)  can  affect  quantitative  predictions 
for  twin  nucleation  [5].  As  shown  later  in  Section  3.4.3,  stability  is  affected  by  the  derivative 
of  ^  at  =  0. 

The  present  variational  framework  for  twinning  and  fracture  does  not  account  for  effects 
of  loading  rate,  temperature  or  explicit  dissipation.  Mathematically,  it  is  sufficient  to  state 
that  the  theory  is  limited  to  isothermal  equilibrium  solutions;  physically,  this  corresponds 
to  neglect  of  local  dynamics  that  presumably  occur  on  very  short  time  scales  for  problems 
of  interest.  Transients,  rate  effects  and  path  dependence  of  solutions  could  be  included  by 
replacing  the  equilibrium  equations  for  order  parameters  by  corresponding  kinetic  equations 
dictating  their  evolution.  For  example,  a  standard  approach  [9,18,20,40,63,64]  involves 
replacing  (2.17)  and  (2.18)  with  Ginzburg-Landau  type  gradient  flow  equations 

90'  90 

W  ”  ^ 

(2.40) 

Here,  the  superposed  dot  denotes  a  material  time  derivative  (X  fixed),  M,,  and  Mj  are 
positive  mobility  coefficients  (inverse  viscosities),  and  0  =  IT  -j-  /  is  total  free  energy 
per  unit  reference  volume.  These  equations  can  be  easily  generalized  to  account  for  kinetic 
coupling  via  replacement  of  scalar  mobilities  with  a  symmetric  2x2  tensor  [8,9].  For 
isothermal  problems,  thermal  effects  can  be  included  via  allowance  of  mobilities,  as  well  as 
surface  energies  and  elastic  coefficients,  to  depend  on  temperature  T.  For  problems  wherein 
temperature  changes  occur,  the  free  energy  function  must  be  extended  to  include  explicit 
dependence  on  temperature  change  (i.e.  specific  heat  and  thermoelastic  coupling),  and  an 
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r/  =  -M„  —  =  -M„  - - V  ■ 
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energy  balance  must  be  introduced  relating  changes  in  temperature,  dissipated  energy  from 
inelastic  work  and  viscosity,  and  possible  heat  conduction  [65].  The  present  assumption 
of  static  mechanical  equilibrium  is  also  strong  and  would  be  inappropriate  for  problems 
of  crack  nucleation  and  extension  driven  by  elastic  waves,  for  example,  dynamic  impact 
involving  spall.  When  inertial  effects  are  important,  the  balance  of  linear  momentum  (2.13) 
must  be  extended  to  include  material  acceleration  as  in  [14,66]. 


3.  Analysis  of  simple  shear 

Considered  in  Section  3  is  application  of  the  general  theory  of  Section  2  towards  a  problem 
of  simple  shear  deformation.  Governing  equations  of  Section  2  are  specialized  towards 
this  deformation  state  in  Section  3.1.  The  general  case  wherein  fracture  and  deformation 
twinning  may  occur,  either  sequentially  or  simultaneously,  in  an  element  of  material  sub¬ 
jected  to  simple  shear  loading  is  studied  in  Section  3.2,  along  with  the  degenerate  case  of 
homogeneous  purely  elastic  deformation.  Considered  in  Section  3.3  is  a  reduction  of  the 
general  treatment  of  Section  3.2  to  the  case  where  fracture/damage  may  occur  but  twinning 
is  prohibited.  Considered  in  Section  3.4  is  a  reduction  of  the  treatment  of  Section  3.2  to  the 
case  where  twinning  may  occur  but  fracture  is  prohibited.  Solutions  are  addressed  for  the 
six  scenarios  shown  in  Figure  1 ;  notation  is  described  in  detail  below. 


3.1.  Governing  equations 

Consider  a  domain  of  initially  homogeneous  material  in  the  Xi  —  Xj  plane.  By  construc¬ 
tion,  during  a  deformation-fracture-twinning  process,  assume  that  fields  u,  f ,  and  q  depend 
only  on  the  X  i  coordinate  and  that  displacement  is  further  restricted  to  only  occur  in  the 
X2  direction,  i.e.  u(Xi)  =  v(X)e2,  where  (Xi,  X2,  X^)  =  (X,  Y,  Z)  and  02  is  a  unit  vector 
tangent  to  Y.  The  length  of  the  body  in  the  X  direction  is  L,  such  that  X  e  [0,  L].  The  width 
in  the  Y  direction  does  not  explicitly  enter  the  problem,  but  the  body  is  considered  “wide” 
enough  to  eliminate  boundary  effects  that  may  lead  to  dependence  of  the  solution  on  Y. 

The  order  parameters  for  fracture,  twinning  and  their  spatial  gradients  become  the  scalar 
fields 

§  =  Vf  =  a|/9A  =  §'(A),  VV$  =  =  d^$/dX^  =  $"(X);  (3.1) 

r]  =  q{X),  Vp  =  dq/dX  =  r]'{X),  XXr]  =  =  d^qldX^  =  r]"(X).  (3.2) 

The  simple  shear  deformation  and  the  deformation  gradient  are,  respectively, 

x  =  A,  y^Y  +  viX),  z  =  Z;  (3.3) 

F  =  1  4- yei  0  ci;  y{X)  —  dv/dX  —  v',  J  —  detF  —  1.  (3.4) 

Twinning  is  restricted  to  occur  in  the  Y  direction  with  habit  plane(s)  normal  to  X,  giving 


■  1  0  O’ 

’o’ 

T 

[FHii)]  = 

yo<p(v)  1  0 

[s]  = 

1 

,  [m]  = 

0 

0  0 1 

0 

0 
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(C) 

l>r- 


0 


Localized  twinning,  zero  stress 
u 


Figure  1 .  (colour  online)  Simple  shearing  problem:  (a)  coupled  twinning  and  fracture,  (b)  no  twinning 
or  fracture,  (c)  stress-free  twinning,  (d)  homogeneous  damage,  (e)  localized  complete  fracture 
and  (f)  inhomogeneous  damage.  Order  parameters  for  twinning  and  fracture  are  rj  and  ^ ;  applied 
dimensionless  shear  displacement  is  v. 


The  elastic  deformation  gradient,  elastic  volume  change,  elastic  deformation  tensor  and  its 
trace  are 


1  +  (/  -  yovf  y  -yo(pO 


1  0  0 

y  -  yovir])  1  0 

0  0  1 


detF^  =  y/7''  =  1; 


[C]  = 


y  -  yov 

0 


1  0 
0  1 


trC  =  3  4-  (y  -  yo(p)  ■ 


(3.6) 

(3.7) 


Assume  ^  —  0  such  that  shear  modulus  degradation  function  0(f)  =  (1  — f)^.  The 
twin  boundary  energy  degradation  function  t  obeys  (2.29)  but  is  otherwise  left  unspecified 
at  present.  Strain  energy  density  of  (2.20)  becomes,  for  the  present  simple  shear  geometry 
with  (f  —  (pip), 

W[y(X),  f  (X),  ^(X)]  =  i/r(trC  -  3)  =  i/rod  -  ?)"(y  -  yo<pf-  (3.8) 


It  is  assumed  for  simplicity  that  the  regularization  length  parameter  h  —  I  is  the  same  for 
twin  boundary  thickness  or  fracture  interfacial  thickness.  Then  (2.30)-(2.32)  become 


^o[?(2f)]  =  tw  :  Vf  0  Vf  =  Tl[^\X)f. 


(3.9) 


Downloaded  by  [Army  Research  Laboratory  ADBV]  at  14:06  01  September  2015 


Philosophical  Magazine 


2673 


A  reduction  of  the  fully  three-dimensional  theory  of  Section  2.3  to  the  one-dimensional 
scenario  considered  here  can  be  achieved  by  setting  /J  — >•  oo  in  (2.31)  with  M  =  ei,  such 
that  fracture  energy  contribution  ( V^)  — >•  oo  for  cracks  of  orientations  differ  from  those 
parallel  to  the  Y  axis.  Thus,  potential  crack  geometries  are  restricted  to  those  shown  in 
Figure  1  since  cracks  of  other  orientations  would  tend  towards  infinite  energy.  Aspects  of 
the  twinning  model  in  (2.26)-(2.28)  reduce  to 

Mri(X),  t(X)]  =  (12r/l)ri\l  -  rifi,  at  :  Vt?  ®  Vt?  =  (3ri/4)[,j'{X)fl  (3.10) 

The  total  energy  functional  (2.9)  becomes,  per  unit  cross-sectional  area, 

'Y  —  j  (W  +  fo  +  go  +  K  :  Xr]  ^  Xr]  +  (o  :  ^X^)dX 

Jo 

i/xo(l  -  -  yovf  +  (12r//)tp2(i  _  ^)2  +  ^^2/^  ^  lrriW{X)f 

+  Tl[^\X)fYx.  (3.11) 

From  (2.13)  and  (2.20)  with  7  =  1,  the  stress  tensor  is  P  =  P-{F^  —  .  It 

follows  that  the  only  nonzero  components  of  the  first  Piola-Kirchhoff  stress  are 

P(X)  =  Pj,Y  =  Pyx  =  fj.oii-^fiy-yo(p),  Pvy{X)  =  -p-o(i-^f(y-yo(p)yo(p- 

(3.12) 

The  equation  for  stress  equilibrium  in  (2.13)  (V  ■  P  =  0)  results  in  the  single  nontrivial 
equation 


P'(X)  =  9P/aX  =  0  =J^  P  =  constant.  (3.13) 

The  conjugate  driving  force  for  fracture  introduced  in  the  first  of  (2.34)  is 

g{X)  =  dw/d^  =  -/rod  -  f)(y  -  yovf.  (3.14) 

The  Euler-Lagrange  equation  for  fracture  order  parameter  equilibrium,  also  in  (2.34), 
reduces  to 

Mo(l-f)(y-yo^)"-2T^//  +  2Tir'-<'[3n(p')V4+12rp2(i_^)2//]  =0.  (3.15) 

The  conjugate  driving  force  for  twinning  introduced  in  the  first  of  (2.33)  is 

r(X)  =  dW/dr]  =  -/royo(l  -  ^fiy  -  yovW ,  (3.16) 

wherecp'  =  d^/d??. Twinning  orderparameterequilibrium,i.e.  the  second  of  (2.33),  requires 

Foyo(i  -  Hf{y  -  yovW  +  Imih”  +  -  24r(p(i  -  3/7  +  2p2)//  =  0.  (3.17) 

Finding  simultaneous  solutions  to  (3.13),  (3.15),  and  (3.17),  given  yet  to  be  prescribed 
boundary  conditions,  becomes  more  straightforward  upon  transformation  of  variables  to 
normalized  dimensionless  form.  Denoting  transformed  variables  by  an  overbar,  define 
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X^XfL,  l^l/L,  V  =  v^fj^o/iTL),  Mo=l,  4' = 'I'L/(2T). 

(3.18) 

These  definitions  result  in 

y  =  dv/dX  =  y^fioL/T,  yo  =  WmqT/T,  P  =  P^L/ipoT)-,  (3.19) 

I' =  =  ^'L,  I"  =  (320) 

i)'  =  dr]/dX  =  r]'L,  i?"  =  d'^rjIdX^  =  ?7"L^.  (3.21) 

The  order  parameters  ^{X)  ^X)  and  ri(X)  — >•  ri(X)  are  unchanged  by  the  transforma¬ 

tion  at  any  material  point  X  =  X/L.  The  elastic  strain  energy  density  becomes 

IT  =  1TL/(2T)  =  i(l  -  t)2(y  -  yocpf.  (3.22) 

Define  the  ratio  R  of  twin  boundary  surface  energy  to  fracture  surface  energy: 

R  =  r/T.  (3.23) 

In  dimensionless  form,  the  total  energy  functional  (3.11)  becomes 

'T(i;,  p,t)  =  -  Yt)V? +rR[6ri^{\  -  rj)'^ /I  +  llW?]  +  n 

4-  {li^'^yx.  (3.24) 

The  Euler-Lagrange  equation  for  stress  equilibrium  (3.13)  becomes 

P'{X)  —  dP/dX  —  0  ^  P  —  (1  —  §)^(y  —  yo(p)  —  constant.  (3.25) 
The  Euler-Lagrange  equation  for  fracture  order  parameter  equilibrium,  (3.15),  reduces  to 
3(1  -  mv  -  vovf  -  H/i-Pil”  -  rRilUp'?  +  6»72(1  -  =  o.  (3.26) 

The  Euler-Lagrange  equation  for  twinning  order  parameter  equilibrium,  (3.17),  becomes 

2(1  -  ?)"yo(y  -  nvW  +  iRm'  +  -  iTiRy^H  -2,n  +  2ri'^)/i  =  0.  (3.27) 

When  the  interpolation  function  (pip)  is  prescribed,  solutions  depend  on  no  more  than  three 
dimensionless  material  parameters:  length  ratio  I,  surface  energy  ratio  R  and  twinning 
eigen-shear  yo- 


3.2.  General  eoupled  problem 
3.2.1.  Boundary  conditions 

Both  ^(X)  e  [0,  1]  and  rjiX)  e  [0,  1]  can  generally  be  nonzero  and  may  possibly  vary 
over  domain  ^2.  Let  v  be  the  applied  shear  displacement  to  the  right  boundary  of  the 
domain  X  e  [0,  1],  with  the  left  end  fixed  with  regards  to  displacement;  i.e.  considered  are 
mechanical  boundary  conditions  of  the  form 

u(l)  =  u  =  f  ydX, 

Jo 


v{0)  —  0, 


(3.28) 
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along  with  vanishing  order  parameter  gradients  at  the  boundaries  [i.e.  r  =  i  =  0  in  (2.15) 
and  (2.19)]: 

|'(0)  =  I'd)  =  0,  i)'(O)  =  i)'(l)  =  0.  (3.29) 

Several  candidate  solutions  of  (3.25),  (3.26)  and  (3.27)  are  now  analysed  for  these  boundary 
conditions. 


3.2.2.  Homogeneous  elastic  deformation 

Perhaps  the  simplest  candidate  solution  is  homogeneous  elasticity,  wherein 

p  =  §  =  0,  P  =  K,  vp(y,o,0)  =  'LE=  5^2.  (3.30) 

Stress  equilibrium  (3.25)  results  in  y  =  u  =  constant.  Twinning  order  parameter  equilib¬ 
rium  (3.27)  at  nonzero  y  requires  either  yo  =  0  or  ^'(0)  =  0,  the  latter  of  which  is  satisfied 
by  choice  (2.35)  in  (2.37),  for  example.  Fracture  order  parameter  equilibrium  (3.26)  requires 
I  =  jdl  —  |)y^  =  0,  which  holds  only  at  null  applied  deformation  (O  =  y  =  0)  or  for 
I  —>■  0.  This  is  the  same  conclusion  reported  later  in  (3.65). 

3.2.3.  Homogeneous  damage  followed  by  localized  fracture  or  localized  twinning 

Next  consider  homogeneous  damage  followed  by  localization  to  either  fracture  (Figure  1  (f)) 
or  twinning  (Figure  1(a)).  Prior  to  localization,  no  twinning  occurs  and  f  is  spatially  constant: 

^{X)  =  =  constant,  r]{X)  —  0.  (3.31) 

Combining  with  equilibrium  leads  to  constant  shear  stress  P  proportional  to  constant  shear 
strain  y  =  F: 

P  =  (1  —  —  constant.  (3.32) 

Twinning  order  parameter  equilibrium  (3.27)  reduces  to  (1  —  yoY (pfO)  —  0  ,  which  at 
nonzero  y  and  <  1  requires  either  yo  =  0  or  ^'(0)  =  0,  the  latter  satisfied  by  choice 
(2.35)  in  (2.37).  Prior  to  localization,  fracture  order  parameter  equilibrium  (3.26)  requires 

(1  -  §h)k"  -  2?h//  =  0,  (3.33) 

which  is  identical  to  (3.61)  of  Section  3.3.2,  with  algebraic  solution 

=  y^/{2/I  +  f),  P  =  y/(l  +  If/2f  (3.34) 

and  total  energy 

=  +  (3.35) 

Therefore,  it  follows  that  the  homogeneously  damaged  state  in  Figure  1(d)  is  a  feasible 
equilibrium  solution  consistent  with  the  governing  equations  of  the  fully  coupled  theory.  A 
lower  bound  for  stability  for  this  homogeneous  solution  is  derived  later  in  Section  3.3.3, 
with  instability  possible  if  y  >  y/l/ifl).  However,  transition  to  a  localized  stress-free  and 
fully  fractured  state  (Figure  1(e))  becomes  energetically  favorable  if  'Fh  >  *Ff  1,  which 
occurs  for  y  ^  2  for  small  I  as  shown  later  in  Figure  3(b).  That  such  a  locally  fully  failed 
state  corresponds  to  a  viable  solution  of  the  coupled  boundary  value  problem  is  verified 
later  in  Section  3.3.4.  However,  transition  from  a  homogeneously  damaged  condition  to 
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this  state  may  require  reversibility  of  damage  in  regions  outside  the  localized  damage  zone, 
which  is  physically  unrealistic. 

Upon  transition  from  homogeneous  damage  to  localized  twinning.  Equations  (3.25)- 
(3.27)  become 


P  =  (1  —  §H)^(y  —  YoV)  —  constant,  (3.36) 

2(1  -  §H)(y  -  -  i'mRilUri'?  +  =  0,  (3.37) 

5(1  -  HnfvoiY  -  YOVW  +  iRlKHtiW  -  12P^(1  -  3r)  +  2ri'^mn)/i  =  0.  (3.38) 


Simultaneous  solution  of  the  above  three  differential  equations  for  y  (2f)  with 

nonzero  shear  stress  P  requires  advanced  numerical  methods  beyond  the  present  scope. 
However,  for  the  scenario  shown  in  Figure  1(a)  with  vanishing  shear  stress  {P  —  0),  (3.36) 
is  trivially  satisfied,  and  (3.37)  and  (3.38)  reduce  to 


-  IhA -  ri^H) R[lhri')^  +  -  7)^/1 


=  0, 


lRmH)ii"]  -  l2Rl^(l  -3>^  +  2n^)i(^^)/l  =  0. 


(3.39) 

(3.40) 


Solution  of  (3.39)  is  difficult  unless  ^  0  and  ('(^h)  =  0,  which  is  the  case  studied 
later  in  Section  3.4.4.  However,  a  stress-free  solution  satisfying  the  remaining  governing 
equations  can  still  be  derived  as  follows.  For  R  >  0  and  <  1,  (3.40)  can  be  reduced  to 
(3.135): 

7j"  -  (32/F)ri^  +  i4S/F)r]^  -  (16/P)r]  =  0,  (3.41) 

which  has  the  inverse  solution  derived  and  evaluated  in  Section  3.4.4: 

X(r^)  =  I  f  dr,/[4ri(l  -  ri)].  (3.42) 

Jo 

The  total  energy  of  the  stress-free,  locally  twinned  domain  becomes,  using  (3.140)  in 
(3.24), 

'I'L  =  £  [t(^u)R[6riFl  -  r)f/l  +  +  ^t^AjdX  =  L(^Yi)R  +  (3.43) 

where  the  first  term  following  the  final  equality  accounts  for  twin  boundary  energy 
degraded  by  damage,  and  the  second  accounts  for  homogeneous  damage.  In  the  absence  of 
homogeneous  damage,  i  —  \  and  total  energy  (which  then  equals  twice  the  normalized  twin 
boundary  surface  energy)  now  takes  the  value  of  R  rather  than  unity  of  (3. 140)  because  of 
the  difference  in  normalization  factors  used  in  (3.105)  and  (3.18).  The  homogeneously 
damaged,  stressed  and  twin-free  domain  and  the  homogeneously  damaged,  unstressed 
domain  with  localized  twinning  have  equal  total  energy  when  (3.35)  and  (3.43)  are  the 
same: 

4'h  =  'Ll  44  1(1  -  ?H)2y2  =  Rii^n).  (3.44) 

Note  that  if  ((§)  =  0(|)  =  (1  —  ^)^,  this  equality  is  independent  of  the  value  of  i-e. 
it  then  does  not  depend  on  the  existence  of  homogeneous  damage  in  the  material.  For  this 
choice  of  i,  the  ratio  of  energy  for  homogeneous  damage,  'I'h  of  (3.35),  to  the  energy  for 
localized  twinning,  'Ll  of  (3.43),  is 


'Lh  _  (l-fe)^y"  +  2gg/r  y" 

'Ll  4(l-tH)2p+2t2/r  2/r-t-y2' 


(3.45) 
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i)  0 

Figure  2.  Ratio  of  energy  for  homogeneous  damage  4'^  to  localized  stress-free  twinning  energy  'I'l 
for  (a)  equal  twin  boundary  and  fracture  surface  energies  (R  =  r/T  =  l)  and  various  dimensionless 
characteristic  lengths  /  and  (b)  various  ratios  of  twin  boundary  to  fracture  surface  energies  {R  —  F/T) 
and  fixed  dimensionless  characteristic  length  /  =  0.01.  Dimensionless  shear  displacement  is  u. 
Transition  from  homogeneous  damage  to  localized  stress-free  twinning  is  energetically  favourable  if 
'I'H/'l'L  >  1. 


This  ratio  is  shown  in  Figure  2(a)  for  various  length  parameters  /  and  R  —  and 
in  Figure  2(b)  for  various  surface  energy  ratios  R  and  fixed  I  —  0.01.  A  transition  from 
homogeneous  damage  to  localized  twinning  occurs  when  applied  displacement  v  —  y 
is  such  that  the  ratio  in  the  first  equality  of  (3.45)  exceeds  unity,  which  corresponds  to 
intersection  of  the  curves  in  Figure  2(a)  or  (b)  with  the  dotted  line  denoted  “twin  threshold”. 
The  intersection  point  is  insensitive  to  I  as  evident  from  Figure  2(a),  but  logically  occurs 
at  lower  applied  displacement  v  for  smaller  R,  which  corresponds  to  relatively  lower  twin 
boundary  surface  energy  (Figure  2(b)).  For  the  choice  t(f )  =  0(f)  used  here,  or  any  other 
choice  with  fn  =  0  and  t(0)  =  1,  transition  from  homogeneous  deformation  to  localized 
stress-free  twinning  becomes  energetically  favourable  for 

u  —  y  >  2\/~R.  (3.46) 

This  can  be  compared  with  the  criterion  derived  later  in  Section  3.3.4  for  transition  from 
homogeneous  elastic  deformation  to  localized  stress-free  fracture: 

0  =  y>2  (I<  0.01).  (3.47) 

Comparison  of  the  two  criteria  leads  to  the  conclusion  that  localized  twinning  should 
occur  at  earlier  applied  shear  displacement  D  than  localized  mode  II  fracture  for  R  <  1, 
and  the  converse  for  R  >  1,  if  the  state  with  lowest  total  energy  is  preferred.  It  should  be 
recalled,  however,  that  the  state  with  total  energy  (3.43)  is  generally  a  nonequilibrium  state 
with  regards  to  the  fracture  parameter  f  since  (3.39)  is  not  necessarily  always  satisfied  for 
the  fully  coupled  model  when  ffj  >  0.  However,  recalling  that  (3.39)  is  satisfied  if  fn  =  0 
[so  long  as  ('(0)  =  0],  the  above  transition  criteria  are  fully  valid  for  evaluating  abrupt 
localized  fracture  or  localized  twinning  from  an  initially  undamaged  state. 
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3.3.  Fracture  in  simple  shear 

Considered  in  Section  3.3  is  reduction  of  the  general  theory  of  Section  3. 1  to  the  case  where 
fracture/damage  may  occur  but  twinning  is  prohibited.  In  terms  of  order  parameters,  this 
implies  r](X)  =  OVZ  e  while  ^(X)  e  [0,  1]  can  vary  over  domain  f2.  In  terms  of 
material  properties,  twinning  would  be  prohibited  by  setting  yo  0  and  F  very  large,  such 
that  4'  of  (2.9)  would  be  penalized  unconditionally  for  any  nonzero  rj  field. 


3.3.1.  Reduced  governing  equations 

Since  no  twinning  occurs,  ^(p)  =  ^(0)  =  0,  F^{r\)  —  F’^0)  —  1  and  F  —  F^  everywhere 
in  Thus,  the  elastic  deformation  tensor  of  (2.6)  and  (3.7)  is,  in  matrix  form,  and  with 
trace:  _ 

1  +  y  0 


[C]  = 


y 

0 


1  0 
0  1 


trC  =  3  +  yL 


(3.48) 


As  in  Section  3.1,  assume  ^  —  0  such  that  0(|)  =  (1  —  |)^.  Strain  energy  density  of 
(2.20)  or  (3.8)  then  reduces  to 


W[Y{X),  t(X)]  =  i/r(trC  -  3)  =  ^/xod  -  ^ 


(3.49) 


Again  denoting  the  length  parameter  h  by  I,  (2.30)-(2.32)  combine  and  reduce  to  (3.9). 
The  energy  functional  (2.9)  or  (3.11)  becomes,  per  unit  cross-sectional  area, 

4'(u,§)  =  j^  c[i/xo(l-?)V"  +  nV^  +  T/(r)"]dX.  (3.50) 

From  (3.12),  or  (2.13)  and  (2.20)  with  7  =  1,  stress  obeys  P  —  ii{F  —  F~'^),  with 
nonzero  components 

P(X)  =  P^Y  =  PyX  =  M0(1  -  ^fy.  (3.51) 

The  Euler-Lagrange  equation  for  stress  equilibrium,  (3.13),  applies  where  now 
specifically 

P'{X)  =  dP/dX  =  0  =►  P/jJLo  =  (1  -  H  fy  =  constant.  (3.52) 


The  conjugate  driving  force  for  fracture  of  (2.34)  or  (3.14)  is 

^(X)  =  aw/9^  =  -Mo(l-t)y".  (3.53) 


The  Euler-Lagrange  equation  for  order  parameter  equilibrium,  (2. 34)  or  (3 . 15),  reduces 
to 

/xo(l-§)y^-2T^/l  +  2Tl^"  =  0.  (3.54) 

Similar  to  [14,31],  apply  definitions  in  (3.18)-(3.20).  The  total  energy  in  (3.24)  or  (3.50) 
becomes  ^ 

'i'(u,|)  =  j|  [i(l-§)2y2+  1^2/;-^  (3  55) 

The  Euler-Lagrange  equation  for  stress  equilibrium  (3.25)  or  (3.52)  becomes 

p'iX)  =  aP/9X  =  0  =J>  P  =  (1  -  ^)2y  =  constant.  (3.56) 
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The  Euler-Lagrange  equation  for  order  parameter  equilibrium  (3.26)  or  (3.54)  likewise 
becomes 

i(l-§)y2_§/[+[|"  =  0.  (3.57) 

Solutions  now  depend  on  only  one  material  parameter,  the  dimensionless  length  ratio  1. 


3.3.2.  Homogeneous  elasticity  and  damage 

Solutions  to  the  problem  outlined  in  Section  3.3.1  are  now  sought  for  which  damage  is 
spatially  homogeneous: 

^(X)  =  =  constant,  |'(X)  =  |"(X)  =  0.  (3.58) 

A  subset  of  this  class  of  solutions  is  shown  in  Figure  1(b),  for  which  the  domain 
remains  perfectly  elastic  and  §  =  OVX  e  [0,  1];  an  example  of  homogeneous  damage  with 
0  <  ^  <  1  is  shown  in  Figure  1(d).  Stress  is  constant  by  (3.56);  it  follows  from  (3.19)  and 
(3.58)  that  strain  y  is  also  constant: 

P  =  (1  —  §h)^K  =  constant  y  =  P/(l  —  §h)^  =  constant.  (3.59) 

Let  the  domain  be  fixed  at  X  =  0,  and  denote  by  v  the  shear  displacement  at  X 

=  Y  [  dX  =  y. 

Jo 

Using  (3.58),  order  parameter  equilibrium  (3.57)  requires 

(1  -  =  2|h//, 


i3(0)  =  0, 


=  F  =  / 


u(l)  —  V  —  j  ydX 


=  1: 

(3.60) 

(3.61) 


which  has  the  algebraic  solution 

=  y^/(2/r  +  y2),  p  =  y/(l  +  ry^/2f.  (3.62) 


Dimensionless  stress  P  is  shown  in  Figure  3(a),  noting  from  (3.60)  that  F  =  y  in  this 
case.  Under  load  control  (i.e.  imposed  nonzero  P),  the  solution  is  nonunique  except  at  the 
load  corresponding  to  maximum  stress  since  two  strain  states  can  produce  the  same  stress 
state.  Stress  attains  a  maximum  value  Pc  at  shear  yc- 


dP/dy  =  0  ^  yc  =  ^2/(30,  Pc  =  ^7^' 


(3.63) 


The  maximum  stress  is  achieved  at  |  regardless  of  /,  so  long  as  I  >0.  Total 

energy  (3.55)  becomes 


'Lh  =  'L[y,^H(K)]  = 


dX=  i(l-§H)"K"+^§f)A.  (3.64) 


2,-,2  ,  lt2  , 


At  maximum  stress  conditions  (3.63),  T'h  =  1/(80-  Total  energy  'I'h  is  shown  in 
Figure  3(b);  as  will  be  derived  later  in  Section  3.3.4,  the  total  energy  associated  with  localized 
stress-free  fracture-see  Figure  1(e)  and  the  dotted  horizontal  line  in  Figure  3(b)  -  is  'Ll  ^  1 
for  I  <  0.1,  regardless  of  applied  displacement.  As  the  characteristic  length  shrinks  to  zero, 
the  purely  elastic  homogeneous  solution  is  recovered: 


(tH,F,'T)^  (0,y,yV4)  as  /  ^  0. 


(3.65) 
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Figure  3.  Homogeneous  analytical  solution  for  (a)  dimensionless  stress  P  vs.  dimensionless 
shear  displacement  u  and  (b)  dimensionless  energy  iFh  vs.  D.  Dimensionless  characteristic 
length  (interfacial  width)  is  Locajized  stress-free  fracture  becomes  energetically  favourable  to 
homogeneous  damage  when  iFh  >  iFp  ^  1  • 


3.3.3.  Stability 

Stability,  or  lack  thereof,  of  the  homogeneous  solution  derived  in  Section  3.3.2,  is  now 
analysed.  The  present  analysis  extends  previous  studies  [14,30,32-35],  and  most  closely 
resembles  the  dimensionless  analysis  in  [31].  In  contrast  to  these  prior  works  that  focused 
on  stability  of  models  involving  linear  elasticity  and  one-dimensional  tensile  loading,  the 
present  analysis  considers  nonlinear  elasticity  and  finite  simple  shear. 

A  family  of  candidate  solutions  is  now  introduced,  where  a  is  a  scalar  parameter  (i.e. 
perturbation)  whose  magnitude  denotes  deviation  from  the  homogeneous  solution  fn- 

^(X,0)  =  §h;  v  =  v{X,a),  y  (X,  a)  =  du/dX.  (3.66) 

Displacement  boundary  conditions  are  imposed,  and  candidate  solutions  are  required 
to  obey  the  same  such  boundary  conditions  imposed  for  the  homogeneous  solution: 

u(0,  a)  =  u(0,  0)  =  0,  v(l,a)  =  v(l,Q)  =  u=  f  y(X,a)dX.  (3.67) 

Jo 

The  candidate  solutions  for  ^  are  further  expected  to  have  vanishing  order  parameter 
gradients  =  d^(X,  a)/dX  at  the  endpoints  (such  a  requirement  is  trivially  satisfied  by 
the  homogeneous  solution): 

f'(0,a)  =  |'(l,a)  =  0.  (3.68) 


Define  the  following  integral  function  of  parameter  a: 


i7(a)=/’  dXI[\-^{X,a)f. 

Jo 

Recalling  from  (3.51)  and  (3.56)  that  shear  stress  obeys 
independent  of  X, 


f 


u=/  y  (X,  a)dX  ^  P(a)  dX /[\  -  ^{X ,  a)Y  ^ 


/' 


(3.69) 

P  =  (1  —  ^Yy  and  is 
P{a)§{a).  (3.70) 
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The  following  relations  then  can  be  derived  directly: 


P{a)  —  v/-&{a),  y  ,  a)  =  u/{i?(o!)[l 
d#(a)  ^  2 

da  Jo  [l-^(X,a)]3  da 
d"i?(a)  _  2  d^$(X,a)^- 

da2  Jo  [l-|(X,a)]3  9a2 


^{X,a)f}- 

'dl(X,ay 

Jo  [l-§(X,a)]4 

da 

(3.71) 

(3.72) 

dX. 


(3.73) 


The  total  energy  of  (3.55)  in  terms  of  parameter  a  is 
'T(a)  = 


iPY  +  k^/I+lh^'f]dX 


02/(2!?)  +  f/l  +  /'(fO^ldX.  (3.74) 


The  first  variation  of  total  energy  is,  upon  integration  by  parts  and  use  of  boundary 
conditions  in  (3.68), 


3'l'(a)  = - ^  /  — d2f5a 

''  "  da 


02 


4!?2 

’2!?2 


C^^dXda^  C  ii-WY±dX,aUW^±\ 
Jo  da  Jo  V/  J  9«  V  9«J 


r^^-idXSa+f 

Jo  (1  -  §)3  aa  Jo 


§)3  da 

Noting  for  the  homogeneous  solution  (a  =  0)  that 


4  -  II"  )  —dXSa. 


da 


(3.75) 


i"iX,  0)  =  0,  !?(0)  =  1/(1  -  §H)^  (3.76) 

and  consulting  (3.61),  the  first  variation  (3.75)  vanishes  identically,  as  it  should  for  an 
equilibrium  state: 

f  f(i-^^)^dXSa+  [  ^^dXaa  =  0.  (3.77) 

2  Jo  da  Jo  I  da 


Appealing  to  (3.72)  and  (3.73),  the  second  variation  of  (3.74)  at  a  =  0  is  computed  as 


S2q/ 

i5a2 


(0)  =  21)2 


dX.  (3.78) 


As  discussed  in  [30],  ^t(O)  >  0  is  necessary  for  stability  of  the  homogeneous  solution, 

and  (0)  >  0  is  sufficient  for  stability  of  the  homogeneous  solution.  Noting  that  the  first 
two  terms  on  the  right  of  (3.78)  are  always  nonnegative,  and  noting  for  the  homogeneous 
solution  that  b  =  y,  a  lower  bound  for  the  unstable  domain  in  terms  of  applied  shear  strain 
is  determined  by  the  third  term  on  the  right  of  (3.78): 


y  >  72/(3/)  =  yc-  (3.79) 

Accordingly,  the  homogeneous  solution  can  become  unstable  under  Displacement-controlled 
loading  when  maximum  uniform  stress  Pq  of  (3.63)  is  exceeded.  Exact  stability  conditions 
for  the  analogous  tensile  problem  are  derived  in  [35],  where  it  is  shown  that  terms  analogous 
to  the  first  two  integrals  on  the  right  of  (3.78)  become  significant  only  for  relatively  large 
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I  >  0.4,  in  which  case  the  maximum  strain  possible  in  the  stable  domain  can  be  larger. 
Otherwise,  for  the  more  physically  realistic  case  of  thinner  crack  boundaries  (i.e.  for  small 
1),  (3.79)  becomes  exact. 


3.3.4.  Localized  complete  fracture 

Considered  next  is  the  problem  shown  in  Figure  1(e):  a  stress-free  domain  with  localized 
damage  concentrated  at  the  midpoint.  Applied  displacement  v  —  u(l)  may  take  any  value 
since  y  may  be  unbounded  without  energetic  consequence  where  the  material  is  fully 
fractured/degraded.  Mechanical  loading  conditions  are 

P  =  [l-|(X)]2y(X)  =  0,  i;(0)  =  0.  (3.80) 

This  implies  y(X)  =  OV^(X)  <  1.  Boundary  and  symmetry  conditions  on  the  order 
parameter  are 

|'(0)  =  I'd)  =  0,  f(2)  =  l.  (3.81) 

The  vanishing  order  parameter  gradient  boundary  condition  implies  conjugate  traction  s  =  0 
in  (2.15)  and  (2.19).  Stress  equilibrium  is  trivially  satisfied.  Order  parameter  equilibrium 
(3.57)  becomes  the  linear,  homogeneous  second-order  differential  equation 

I"  -  f /F  =  0,  (3.82) 

which  has  the  immediate  solution 

^(X)  =cie^/'  +  C2e-^/';  (3.83) 

Cl  =  C2  =  1  -b  for  X  <  i 

Cl  =  l/[e'/(20  +  e3/(20],  C2  =  +  e-'/^^')]  for  X  >  i.  (3.84) 

The  dimensionless  total  energy  functional  in  (3.55)  is  then 

fi  .  .  /'1/2  _  __  _  e^/'-t-e”'/' 

q/p  =  1  /  =  /  [^/l  +  lih^]dX  =  (3.85) 

Jo  Jo  [ei/(20 +  e-i/(2/)]2 

A  direct  calculation  verifies  that  'Tp  ^  1  for  I  <0.1,  as  depicted  by  the  horizontal  dotted 
line  in  Figure  3(b).  As  is  evident  from  this  figure,  for  small  I,  the  localized  stress-free 
solution  has  greater  total  energy  than  the  homogeneous  damage  solution  ('Fp  >  'Fh)  for 
applied  shear  F  <  2.  However,  stability  of  the  homogeneous  elastic  solution  in  the  context 
of  the  perturbation  analysis  of  Section  3.3.3  holds  for  0  <  yc  —  ^/2/(3l).  For  example, 
taking  I  =  0.01,  stability  holds  for  F  <  8,  but  the  homogeneous  solution  of  Section  3.3.2 
represents  a  local  minimum  energy  configuration  (i.e.  is  in  fact  metastable)  for  2  <  F  <  8. 

Verification  is  straightforward  that  this  candidate  solution  involving  stress-free,  local¬ 
ized  fracture,  in  the  absence  of  twinning,  is  valid  for  the  fully  coupled  theory  of  Section 
3.2.  Since  rj(X)  =  0,  P  =  (1  —  §)^y  —  0  ^  y  —  OV^  <  1.  Twinning  order  parameter 
equilibrium  (3.27)  reduces  to  (1 —§)^yoK<p'(0)  =  0,  which  is  identically  satisfied  by  P  =  0. 
Fracture  order  parameter  equilibrium  (3.26)  reduces  to  (3.82);  the  solution  for  §  (X)  is  (3.83) 
with  integration  constants  (3.84).  The  total  energy  of  this  stress-free  domain  with  localized 
complete  fracture  (^  =  1  at  X  =  j)  is  given  by  (3.85),  i.e.  in  notation  of  Section  3.2.3, 
qfp  =  %  1  for/"  <  0.1. 
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3.3.5.  Inhomogeneous  damage 

Solutions  are  sought  to  the  shear  problem  wherein  damage  (i.e.  § )  may  be  nonuniform  over 
domain  [0,  1];  see  Figure  1(f).  The  solution  procedure  is  similar  to  that  described  elsewhere 
for  tensile  loading  [14,31,32]. 

The  loading  parameter  is  stress  P,  constant  over  the  domain  according  to  (3.56). 
Boundary  conditions  are  imposed  so  that  the  homogeneous  solution  of  Section  3.3.2  is 
also  a  solution  to  the  present  problem: 

§(0)  =  §(1)  =  l'(0)  =  I'd)  =  0.  (3.86) 


Here,  fn  (F)  is  the  value  of  the  order  parameter  in  the  corresponding  homogeneous  solution 
during  the  loading  phase  prior  to  instability  (^h  <  5)  is  obtained  by  simultaneous 
solution  of  (3.62)  with  P  <  Pq  imposed.  The  following  symmetry  conditions  are  also 
imposed: 

^'(5)=®’  ^(5-^)  =1(2+^)-  (3-87) 

The  displacement  at  the  right  end,  i;(l)  =  v,  can  be  determined  as  an  outcome  of  the 
solution. 

The  Euler-Lagrange  equation  for  order  parameter  equilibrium  (3.57)  becomes,  in  terms 
of  P  =  (1  —  §)^y  =  constant  and 

3p2/(l-?)^-§/l+l|"  =  0.  (3.88) 

Its  solution  is  facilitated  by  the  following  change  of  variables: 

5(X)  =  1-§(X),  d(X)  =  -|'(X),  s"{X)  =  -l"(X).  (3.89) 


In  terms  of  s  and  P,  (3.88)  becomes 


d  /  1 


Using  the  identities 


Equation  (3.90)  can  be  written 

d 
di 


Pis'f 

2 


P^J/s 

1  =  0. 

1 

d 

r(^')"i 

2^3’ 

2 

4^2  ) 

+  d.  1 

.2) 

-d)  =  0. 


(3.90) 


(3.91) 


(3.92) 


Defining  =  sjyi  =  1  —  f  (5)  and  noting  that  ^'(j)  =  0  from  (3.87),  integration  of 
(3.92)  over  domain  [^m,  corresponding  to  right  half  X  results  in 

Fis'f/2  -  P^l/  (4^2)  +  s^/2  -  s  =  -P^l/  (44)  +  4/2  -  iM.  (3.93) 

Defining  SH  =  l—|H,andnoting  boundary  conditionss(l)  =  SHands'(l)  =  (ds/dX)(l)  = 
0,  the  following  equation  can  be  solved  numerically  for  sm  for  a  given  P: 


P^l/  ^44^  “■  4/7  +  ■sm  —  p2//  (^44^  ~  ■^h/7  +  ■sh- 


(3.94) 
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Figure  4.  Inhomogeneous  numerical  solution  for  order  parameter  ^  vs.  dimensionless  coordinate 
X.  Dimensionless  characteristic  length  is  /.  Ratio  of  applied  stress  to  critical  stress  at  instability  is 
X  =  P/PC- 


Inversion  of  (3.93)  leads  to  a  differential  equation  for  X{s): 

dX/ds  =  1/s'  =  (/7V2)  [p2;7(44)  +  4/2  -SM  +  P^I/(4s^)  -  s^ /I  +  ; 

(3.95) 

integration  results  in 

X{s)  =  (772)  [^77(44)  +  4/2  -  +  p77(47)  -  7/2  +  sY^'^  As.  (3.96) 


This  can  be  numerically  evaluated  and  then  inverted  to  provide  the  order  parameter  profile 
^{X)  =  1  —  i(X)  for  X  e  [7  l];  the  profile  for  X  <  j  is  then  obtained  by  applying 
symmetry  conditions  (3.87).  Results  are  shown  in  Figure  4  for  /  =  0.01,  where  x  =  P / Pc 
is  the  applied  stress  normalized  by  critical  stress  of  (3.63).  As  the  applied  stress  decreases 
(increases),  fracture  in  the  inhomogeneous  solution  becomes  more  localized  (diffuse).  The 
total  energy  and  right  side  displacement  are  obtained  by  numerically  evaluating  the  integrals 

l\\PY  +  \H'^/i+\Kh^]AX,  D(P)  =  /oydX  =  /J[P/(l-|)2]dX. 

Jo 

(3.97) 

Results  are  compared  with  those  of  the  homogeneous  solution  of  Section  3.2  (i.e.  4'h 
and  Ufj)  in  Table  1.  For  x  <  1,  the  total  energy  and  right  side  shear  displacement  of 
the  inhomogeneous  solution  exceed  that  of  the  homogeneous  solution;  as  x  — 1,  the 
inhomogeneous  solution  degenerates  to  the  homogeneous  solution. 


3.4.  Twinning  in  simple  shear 

Considered  in  Section  3.4  is  reduction  of  the  general  theory  of  Section  3.2  to  the  case  where 
deformation  twinning  may  occur,  but  fracture  is  prohibited.  This  implies  |  (Z)  =  OVZ  e  ^2, 
while  ri{X)  e  [0,  1]  can  vary  spatially  over  domain  f2.  In  terms  of  material  properties. 
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Table  1 .  Comparison  of  total  energy  and  displacement  between  inhomogeneous  and  homogeneous 
solutions  (nonlinear  phase  field  fracture  model)  under  load  control;  is  the  maximum  value  of  the 
order  parameter  in  the  inhomogeneous  solution,  and  fn  is  the  homogeneous  order  parameter. 


X 

0 

0H 

V 

i4t 

0.35 

1.3 

0.7 

1.8 

1.6 

0.87 

0.01 

0.64 

3.1 

2.5 

3.5 

3.3 

0.72 

0.05 

0.85 

5.8 

5.4 

5.1 

4.9 

0.57 

0.11 

0.97 

9.1 

8.8 

6.7 

6.5 

0.40 

0.18 

1 

12.5 

12.5 

8.2 

8.2 

0.25 

0.25 

fracture/damage  would  be  prohibited  by  setting  T  very  large  and  0  =  ?  =  1,  such  that  4' 
of  (2.9)  would  be  penalized  unconditionally  for  any  nonzero  §  field. 


3.4.1.  Reduced  governing  equations 

Twinning  deformation  (2.5)  is  of  the  form  in  (3.5),  and  the  elastic  deformation  tensor  of 
(2.6)  is  identical  to  that  in  (3.7).  Noting  that  in  the  absence  of  damage,  /r  =  jjLQ,  0  =  0  and 
0  =  t  =  1,  strain  energy  density  of  (2.20)  or  (3.8)  reduces  to 


=  2Mo(trC  -  3)  =  5/ro(y  -  Yovf-  (3-98) 

Equations  (2.26)-(2.28)  or  (3.10)  reduce  to 

/o[t7(2f)]  =  (12r//)t?2(l  _  ^)2,  ;  Vr,  0  Vr,  =  (3r//4)[t?'(X)]2.  (3.99) 

The  energy  functional  (2.9)  or  (3.11)  becomes,  per  unit  cross-sectional  area, 

'T(U,  q)  =  £  -  Yovf  +  urq\l  -  qf/l  +  |r/(p')"]dX.  (3.100) 

From  (3.12),  or  from  (2.13)  and  (2.20)  with  7  =  1,  the  stress  tensor  reduces  to  P  = 
.  The  only  resulting  nonzero  components  are 

P(X)  =  =  Pyx  =  /ro(K  -  YoV),  PyviX)  =  -/xo(k  -  Yo<p)yqV-  (3.101) 

The  Euler-Lagrange  equation  for  stress  equilibrium  (3.13)  is  now 

P'(X)  —  dP/dX  —  0  ^  P/mo  —  Y  ~  YoV  =  constant.  (3.102) 

The  conjugate  driving  force  for  twinning  of  (2.33)  or  (3.16)  reduces  to 

r(X)  =  dWIdq  =  -moKo(K  -  YovW ^  (3.103) 

where  =  dtp/dq.  The  Euler-Lagrange  equation  for  order  parameter  equilibrium  in  (2.33) 
or  (3.17)  is  now 


^riq"  -  24r?7(l  -3q  +  2q^)/l  +  fiQYoiY  “  YovW  =  0. 


(3.104) 
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Similar  to  Section  3.1,  dimensionless  variables  denoted  with  overbars  are  used  in  Section 
3.4;  however,  twin  boundary  surface  energy  F,  rather  than  fracture  energy  T,  is  applied  for 
normalization  here  [(3.21)  still  holds]: 


X^X/L,  l^l/L, 

Y  =  dv/dX  =  yV/roi/r, 


V  =  V/ro/CFL),  MO  =  1,  4'  =  'hL/(2r); 

(3.105) 

Yo  =  yo^/tJ^oL/^,  P  =  PyZ7(Mon  =  K  -  YoV- 

(3.106) 


The  total  energy  functional  (3.24)  or  (3.100)  becomes,  in  normalized  form, 

4'(i;,  p)  =  [i(y  -  Y^qif  +  -  ^)2/i'  +  )^]dX.  (3.107) 

The  Euler-Lagrange  equation  for  stress  equilibrium  (3.25)  or  (3.102)  becomes 

P'(X)  =  dP/dX  =  0  P  =  y  —  yo'?’  =  constant.  (3.108) 

The  Euler-Lagrange  equation  for  order  parameter  equilibrium  (3.27)  or  (3. 104)  becomes 

3;Y  -  12,7(1  -  3p  +  2p2)/[  +  1  yo(y  -  yovW  =  0.  (3.109) 

Solutions  now  depend  on  only  two  material  parameters,  dimensionless  length  ratio  I  and 
twinning  eigen-shear  yo,  when  interpolation  function  (p(ri)  is  fully  prescribed. 


3.4.2.  Homogeneous  elasticity  and  twinning 

Solutions  to  the  problem  outlined  in  Section  3.4.1  are  now  sought  for  which  twinning  is 
spatially  homogeneous: 

ri{X)  =  rju  =  constant,  rj'(X)  =  rj"(X)  =  0.  (3.110) 

A  subset  of  this  class  of  solutions  is  shown  in  Eigure  1(b),  for  which  the  domain  remains 
perfectly  elastic  and  p  =  OVX  e  [0,  1].  Shear  stress  is  constant  by  (3.108).  Since  (pirm)  = 
constant,  it  follows  from  (3.106)  and  (3.110)  that  strain  y  is  also  constant: 

P  —  Y  —  YoVivn)  —  constant  ^  y  —  P  +  YoVi^n)  =  constant.  (3.111) 

Let  the  domain  be  fixed  at  X  =  0,  and  denote  by  v  the  shear  displacement  at  X  =  1: 

=  Y  f  dX  =  y.  (3.112) 

Jo 

Using  (3. 110),  order  parameter  equilibrium  (3. 109)  requires 

24,7(1  -  3,7  +  2,7")  =  /yo(y  -  YovW-  (3.113) 

Eurther  analysis  requires  selection  of  interpolation  function  (p(rj)  first  introduced  in  (2.5). 
Choosing  tp  —  3p^  —  2r\^  from  (2.35),  equality  (3.113)  becomes  the  fifth-order  polynomial 
expression 

6y2p5  -  15y27^  +  (9Ko  +  24//)  ,7^  +  (3y  yo  -  36//)p2  _  (3j7  _  12//)^  =  0.  (3.114) 


u(0)  =  0, 


=  0  =  / 


u(l)  —  V  —  !  ydX 
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This  polynomial  has  two  imaginary  roots  and  the  three  real  roots  listed  below: 

1  3a  2b  d 

ha  =0>  'IH  =  1,  =  :::  +  :^  +  (3.115) 

2  Id  3d  oa 

a^y^,  b^U/I,  c^yyo,  (3.116) 

d  ={27fl^  -  54-a^c  +  2^/a^[16b^  +  243a^(b  -  3c)  -  27fl(4b2  -  27c2)]  (3.117) 


Numerical  evaluation  of  the  third  solution  in  (3.115)  over  physically  realistic  ranges  of 
parameters  10"  ^  <  /  <  0.1,  1  <  yo  <  and  0  <  y  <  Yo  produced  no  values  within 
the  constraints  rjn  e  (0,  1).  Thus,  the  only  feasible  homogeneous  solutions  are  ?7h  =  0  and 
?7h  =  1-  (A  recent  computational  phase  field  study  of  spinodal  twinning  [38],  on  the  other 
hand,  did  find  possible  homogeneous  solutions  with  uniform  order  parameter  values  in  the 
range  0  <  ?7h  <  1,  though  governing  and  boundary  conditions  imposed  in  that  work  differ 
from  those  considered  here.)  The  total  energy  functional  (3.107)  becomes 


'k(u,  ph) 


-f:b 


4(K  -  yo<pr  +  6r]^(l 


rinf/l 


dX. 


(3.118) 


For  =  0  or  PH  =  1,  using  (2.37)  and  noting  here  that  y  —  v. 


^E  =  vp(y,0)  =  .i/T  =  'T(y,l)  =  d(y-yo)2,  (3.119) 


where  tpE  is  the  energy  of  the  elastic  solution  (no  twin)  and  'Tx  is  the  energy  of  the  fully 
twinned  homogeneous  solution.  These  two  energies  are  equal  at  y /yo  —  j,  with  the  former 
purely  elastic  solution  in  Figure  1(b)  energetically  favorable  for  applied  shear  y/yo  <  2 
and  the  latter  fully  twinned  solution  favourable  for  y/yo  >  j- 


3.4.3.  Stability 

Stability,  or  lack  thereof,  of  homogeneous  purely  elastic  solution  p  =  pn  =  0  is  now 
analysed.  Interpolation  function  ^(p)  is  left  general  until  the  conclusion  of  the  analysis. 

A  family  of  candidate  solutions  is  now  introduced,  where  a  is  a  scalar  parameter  (i.e. 
perturbation)  whose  magnitude  denotes  deviation  from  the  homogeneous  solution  pn: 

p  =  p(X,a),  p(X,0)  =  pH  =  0;  v  —  v{X,a),  y(X,a)  —  dv/dX.  (3.120) 

Displacement  boundary  conditions  are  imposed,  and  candidate  solutions  are  required  to 
obey  the  same  such  boundary  conditions  imposed  for  the  homogeneous  solution: 

£i(0,  a)  =  u(0,  0)  =  0,  i}(l,  a)  =  u(l,  0)  =  b  =  f  y(X,a)dX.  (3.121) 

Jo 

Candidate  solutions  for  p  are  further  expected  to  have  vanishing  order  parameter  gradients 
p'  =  dri(X,  a)/dX  at  the  endpoints,  as  trivially  satisfied  by  the  homogeneous  solution: 

p'(0,Q!)  =  p'(l,Q!)  =  0.  (3.122) 

Define  the  following  integral  function  of  parameter  a: 

d(a)  —  f  (pdX,  <p  —  (p[ri(X ,  a)],  dd/da—  f  (p'{drj/da)dX.  (3.123) 

Jo  Jo 
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Recalling  from  (3.106)  that  shear  stress  obeys  P  —  y  —  y^cp  and  is  independent  of  X, 

+  yo  f  (pdX  =  P(a)  +  yod(a),  (3.124) 

Jo 

y(X,  a)  =  P(a)  +  yo(p[r]iX,  a)]  =  D  +  yoi(p  -  d).  (3.125) 

The  total  energy  of  (3.107)  in  terms  of  parameter  a  then  follows  as 

«i/(a)  =  i  [  -  2Dyo‘P^  +  yid^)dX  +  3  [  -  r^f/I  +  ll^fjdX.  (3.126) 

Jo  Jo 

The  first  variation  of  total  energy  is,  upon  integration  by  parts  and  use  of  boundary  conditions 
in  (3.122), 


/' 


V  —  I  y(X,a)dX  —  P(a) 


S^){a)  =  - 


m 
a 


__  d#  __  dr]' 

(Ko  ^  ~  ^yov):, - vyodrp  — 

da  da 


dX&a 


1  r4 

,,(i 


3ri^  +  2r]^)  -  -r]” 


dn  - 
— dX&a. 
da 


Noting  for  the  homogeneous  solution  (a  =  0)  that 


(3.127) 


rjiX  ,0)  —  Tj^^X  ,0)  —  0,  d  —  rp  —  O,  d?7/da  =  ^'9?7/9a,  (3.128) 

the  first  variation  of  total  energy  evaluated  at  a  =  0  vanishes,  as  it  should  for  an  equilibrium 
state: 


9'T(0)  =  0. 


(3.129) 


The  second  variation  of  energy  functional  (3.126)  at  a  =  0  is  computed  upon  successive 
differentiation  as 


Noting  that  the  rightmost  term  of  (3.130)  is  always  nonnegative,  and  noting  that  for  the 
homogeneous  solution  that  u  =  y,  a  lower  bound  for  the  unstable  domain  in  terms  of 
applied  shear  strain  is  determined  by  the  sum  of  coefficients  of  the  first  two  integrals  on  the 
right  of  (3.130): 


-  yo  ,  12 

y  >  ~ — f  = - r. 

2  /yo[^'(0)]2 


(3.131) 


If  <p'(0)  =  0,  as  occurs  in  (2.37)  for  particular  choice  (2.35)  considered  in  Section  3.4.2, 
the  homogeneous  elastic  solution  is  stable  with  respect  to  perturbation  a  (i.e.  |;^(0)  >  0) 
for  any  applied  strain  since  I  >  0.  Even  if  stable  in  the  above  sense,  such  an  elastic  solution 
may  reflect  a  local  minimum  energy  configuration  (i.e.  metastability),  since  other  solutions 
(e.g.  a  heterogeneous  rj  field  or  fully  twinned  with  =  1)  may  have  lower  total  energy 
at  the  same  imposed  v.  For  a  linear  interpolator  (e.g.  cp  =  rj)  similar  to  that  considered  in 
[40],  stability  of  the  homogeneous  elastic  solution  is  not  always  ensured  for  y /yo  > 
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3.4.4.  Localized  complete  twinning 

Considered  next  is  the  problem  of  Figure  1  (c):  a  stress-free  domain  with  localized  twinning  at 
its  midpoint.  Applied  displacement  v  —  i;(l)  takes  a  value  within  0  <  0  <  yo-  Mechanical 
loading  conditions  are 

P  =  y(X)  -  yo^[?7(X)]  =  0  y  =  yo?),  u(0)  =  0.  (3.132) 

Boundary  and  symmetry  conditions  on  the  order  parameter  are 

p'(0)  =  p'(l)  =  0,  *7(2)  =  1,  =  p(i+X).  (3.133) 

Vanishing  order  parameter  gradients  at  the  endpoints  corresponds  to  r  =  0  in  (2.15)  and 
(2.19).  Order  parameter  equilibrium  (3.109)  can  be  written  as 

+  \Pyqv'  -  (24/0773  +  (36/0772  -  (12/077  =  0.  (3.134) 

Substituting  the  null  stress  condition  (3.132),  this  becomes  the  nonlinear  second-order 
differential  equation 

7)"  -  (32/P)773  -P  (48/P)772  -  (16/P)77  =  0.  (3.135) 

Notice  that  for  null  shear  stress,  the  particular  choice  of  function  cp  becomes  irrelevant. 
Differential  equation  (3.135)  can  be  written  as 

d  Ibdn  9 

)  (3.136) 

d77  O  d77 

Integrating  both  sides  of  (3. 136)  from  the  left  side  of  the  domain  where  77  (X)  =  77(0)  =  0  to 
any  77  at  point  X  and  taking  the  positive  root  gives  the  first-order  nonlinear  differential 
equation 

drj/dX  =  477(1  -  77)/L  (3.137) 

Inversion  of  this  equation  and  integration  produce 

X(77)  =  [/'^d77/[477(l-77)].  (3.138) 

Jo 

Integral  (3.138)  is  evaluated  numerically  (here  via  trapezoidal  quadrature)  and  then  inverted 
to  give  order  parameter  profile  r]{X)  over  the  left  half  of  the  domain,  i.e.  over  0  <  X  <  ^. 
Points  Xi  and  Xj  bounding  the  region  where  0  <  77  <  1  are  determined  simultaneously  by 
the  condition  77  =  IVX  e  (X2,  j]  and  the  magnitude  of  applied  shear  that  sets  the  width  of 
the  fully  twinned  region: 

/  ydX  +  yo(l  -  2X2).  (3.139) 

Jxi 

The  order  parameter  profile  for  the  right  half  of  the  domain  can  be  constructed  using 
symmetry  conditions  in  (3.133).  Shown  in  Figure  5(a)  is  the  order  parameter  profile  for 
applied  displacement  u  =  yo/2.  The  domain  shown  is  limited  to  the  relatively  small  region 
where  77  varies  with  X.  For  example,  for  I  —  0.01,  Xi  ^  0.22  and  X2  ^  0.24.  The  width 
of  the  twin  boundary  region  X2  —  Xi  increases  with  increasing  1. 
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X  0/Vo 


Figure  5.  Inhomogeneous  numerical  solution  (a)  for  stress-free  twinning  at  applied  dimensionless 
shear  displacement  of  u  =  yo/2:  profile  of  order  parameter  r]  at  the  left  twin  boundary.  Dimensionless 
coordinate  is  X\  dimensionless  characteristic  length  is  /;  normalized  complete  twinning  shear  is  yq. 
Ratio  (b)  of  homogeneous  elastic  energy  W  to  localized  stress-free  twinning  energy  '1'l  vs.  applied 
normalized  shear  v/yq.  Transition  from  homogeneous  elasticity  to  localized  stress-free  twinning  is 
energetically  favourable  if  W /T'l  >  1. 


The  total  energy  of  the  localized  stress-free  solution  is  found  by  application  of  (3.107), 
which  reduces  here  via  (3.132)  (yielding  strain  energy  W  —  Q  from  P  —Q)  and  (3.137)  to 


'Ll  = 


/'1/2 

2  /  -  nf/l  +  ll{Ar^{\  -  ri)/lf}AX  =  (24//') 

Jo 


rj^il  -  rifdX  =  1. 


(3.140) 

The  final  value  'Ll  =  1  is  consistent  with  normalization  (3.105)  giving  'Ll  =  2r /L  (i.e.  the 
exact  surface  energy  of  two  unstressed  twin  boundaries)  and  was  verified  numerically  for 
lo-f’  <  I  <  0.2.  Shown  in  Figure  5(b)  is  total  energy  for  the  homogeneous  elastic  solution 
of  (3.119)  of  Section  3.4.2,  W  —  'Ll  =  normalized  by  'Ll  =  1.  Localized  stress- 
free  twinning  under  displacement  control  (Figure  1(c))  becomes  energetically  favourable 
to  homogeneous  elasticity  (Figure  1(b))  when  u  is  large  enough  such  that  W  >  1  'Ll  > 
'Ll  44  y  >  2.  The  larger  the  value  of  yo,  the  greater  elastic  energy  relieved  by  twinning  at 
a  given  value  of  u/yo- 

Solution  of  a  more  complex  boundary  value  problem  involving  nonzero  shear  stress  P 
and  heterogeneous  or  localized  twinning  corresponds  to  solution  of  (3. 134)  with  appropriate 
boundary  conditions  and  requires  choice  of  interpolation  function  cp.  Solution  of  this  more 
difficult  differential  equation  should  be  possible,  but  requires  more  advanced  numerical 
methods  beyond  the  present  scope. 

Verification  is  straightforward  that  this  candidate  solution  involving  stress-free,  local¬ 
ized  twinning,  in  the  absence  of  fracture,  is  valid  for  the  fully  coupled  theory  of  Section  3.2. 
Since  §(X)  =0,  P  =  y  —  yo  =  0=J>y  =  yo*?-  Fracture  order  parameter  equilibrium 
(3.26)  requires 


+  6??2(1  -  =  0.  (3.141) 


This  equation  can  be  satisfied  unconditionally  by  choosing  a  degradation  function  satisfying 
('(0)  =  0;  for  example,  a  possible  choice  also  obeying  requirements  in  (2.29)  is  t(f)  = 
1  —  q){^)  =  1  —  rf'Q  —  2rf).  An  even  simpler  example  is  t  =  1,  which  decouples  energies  of 
twin  boundaries  and  fracture  surfaces.  Twinning  order  parameter  equilibrium  (3 .27)  requires 
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((0)R[|Zi7"  -  12)7(1  -  3)7  +  2r]^)/I]  =  0,  (3.142) 

which  is  identical  to  (3.135)  for  finite  R  =  F/T  since  )(0)  =  1  by  (2.29).  Therefore, 
it  follows  that  this  solution,  depicted  in  Figure  1(c),  is  admissible  in  the  fully  coupled 
theory  so  long  as  l  is  chosen  appropriately.  Choosing  t(f)  =  0(^)  is  incompatible  with 
such  a  solution  since  then  T  =  2(f  —  1)  and  (3.141)  is  no  longer  satisfied,  in  general,  for 
nonzero  R. 


4.  Model  predictions  for  real  crystalline  solids 

Upon  consideration  of  stabilities  and  total  energies  of  various  solutions  derived  in  Section  3 
for  the  phase  field  constitutive  model  of  Section  2  applied  to  a  finite  simple  shear  problem, 
various  criteria  for  shear-induced  fracture  and  twinning  can  be  deduced.  These  criteria  can 
be  applied  to  real  materials  by  converting  normalized  quantities  to  their  dimensional  forms, 
as  demonstrated  next. 

First  consider  mode  II  fracture.  The  shear  stress  at  which  a  homogeneously  damaged 
solid  becomes  unstable  according  to  (3.63)  and  (3.79)  is  converted  to  dimensional  form  via 
(3.19),  leading  to 

Pc  =  ^V6moT//.  (4.1) 

The  shear  stress  at  which  localized  fracture  becomes  energetically  favourable  to  homoge¬ 
neous  elastic  deformation  is  found  according  to  derivations  in  Section  3.3.4  and  listed  in 
(3.47): 

Pp  =  lyjfioTIL  for  /  <  O.OIL.  (4.2) 

These  two  criteria  for  shear  fracture  can  be  compared  with  the  so-called  theoretical  fracture 
strength  [3,44,46]  needed  to  relatively  displace  (here  under  shear  loading)  two  slabs  of  linear 
elastic  material  here  by  a  slip  distance  equal  to  the  spacing  between  interatomic  planes: 

Pthjac^  Po/i^TT).  (4.3) 

Unlike  theoretical  strength  (4.3),  which  depends  only  on  shear  modulus,  the  new  criteria 
presented  in  (4.1)  and  (4.2)  are  proportional  to  the  square  root  of  the  product  of  shear 
modulus  and  surface  energy. 

Distinctions  and  similarities  among  the  fracture  criteria  discussed  above  and  Griffith’s 
criterion  [67]  are  noteworthy.  Under  plane  strain  conditions,  in  the  context  of  pure  mode 
II  loading,  Griffith’s  criterion  for  extension  of  a  straight  crack  of  length  2a  in  an  infinite 
isotropic  linear  elastic  medium  can  be  stated  as  [43] 

Pc  =  2y/roT/[(l  -  v)jTa],  (4.4) 

with  Pq  the  far-held  shear  stress  and  v  Poisson’s  ratio.  Griffith’s  criterion  (4.4)  applies  for 
a  pre-existing  sharp  crack  of  hnite  length  in  an  inhnite  body,  while  (4.1)-(4.3)  correspond 
to  nucleation  of  a  crack  in  an  initially  homogeneous  crystal  of  hnite  size.  While  (4.1),  (4.2) 
and  (4.4)  all  incorporate  a  shear  modulus  and  surface  energy,  only  the  latter  (Griffith)  also 
involves  compressibility  via  Poisson’s  ratio.  Note  also  that  (4.2)  and  (4.4)  give  identical 
predictions  when  domain  size  L  =  (1  —  v)7Ta  ^  2a.  Griffith’s  criterion  can  be  extended 
for  applications  to  bodies  of  hnite  size  by  accounting  for  dependence  of  the  crack  tip  stress 
intensity  factor  on  geometry  of  the  body,  in  which  case  the  critical  far-held  stress  for  crack 
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extension  depends  on  both  the  initial  crack  length  and  the  position  of  the  crack  relative  to 
external  surfaces  of  the  body. 

It  is  emphasized  that  all  fracture  criteria  considered  thus  far  apply  for  brittle  elastic  solids. 
When  crack  tip  plasticity  becomes  important,  an  approximate  and  simple  way  of  considering 
the  corresponding  increase  in  toughness  involves  replacement  of  surface  energy  T  in  the 
above  criteria  with  an  effective  (increased)  surface  energy  accounting  for  plastic  work  in 
the  process  zone  [68].  A  more  physically  detailed  and  complete  model  for  ductile  crystals 
would  involve  including  additional  order  parameters  (and  their  corresponding  effects  on 
kinematics  and  energy)  in  the  phase  field  framework  accounting  for  slip  of  full  and/or  partial 
dislocations  on  each  relevant  slip  system  [37,69]. 

Next,  consider  deformation  twinning.  The  shear  stress  at  which  localized  twinning 
becomes  energetically  favourable  to  homogeneous  elastic  deformation  is  found  according 
to  derivations  in  Section  3.4.4  and  listed  in  (3.46): 

Pr  =  2V^y/roT/L  =  ly/pL^T/L.  (4.5) 

This  new  criterion,  like  those  above  for  fracture,  is  proportional  to  the  square  root  of 
the  product  of  a  shear  modulus  and  a  surface  energy.  Criterion  (4.5)  can  be  compared  with 
two  others  from  the  literature,  the  theoretical  strength  associated  with  twinning  [70,71]: 

^’th.twin  =  ftoyo/(2:r).  (4.6) 

and  the  twinning  stress  associated  with  moving  a  partial  twinning  dislocation  of  Burgers 
magnitude  bp  against  the  resistance  associated  with  stacking  fault  energy  Wsp  ^  2r  [3,46]: 

-PsF,  twin  =  2r/bp.  (4.7) 

Representative  properties  for  a  generic  material  of  the  sort  considered  most  often  in 

derivations  of  Section  3  as  well  as  for  three  real  materials-sapphire,  calcite  and  magnesium- 
are  listed  in  Table  2.  Properties  of  the  real  crystalline  solids  are  obtained  from  [5,8,39]  and 
references  quoted  therein,  with  the  exception  of  cohesive  surface  energy  for  magnesium 
[72]  and  bp,  the  latter  which  is  necessarily  obtained  from  additional  references  for  sapphire 
[50],  calcite  [4]  and  magnesium  [73]. 

Each  of  these  crystals  exhibits  multiple  twin  systems  in  nature.  Considered  here  are 
twinning  on  the  rhombohedral  (R)  and  basal  (B)  systems  in  sapphire,  respectively  (1  101) 
(1  102}  and  (1  1  00)[000  1}  in  hexagonal  Miller  indices.  The  e"*"  twin  system  in  calcite 
corresponds  to  (00  1)(1  10}  in  the  cleavage  rhomobohedral pseudocell  [39,42].  The  primary 
tensile  twin  system  in  magnesium  corresponds  to  (101  1)[1012}  in  hexagonal  Miller 
indices  [8]. 

Criteria  (4.1),  (4.2)  and  (4.5)  depend  on  a  length  parameter  (/  or  L),  with  strength 
increasing  as  size  decreases  in  an  inverse  square-root  manner,  a  phenomenon  widely 
observed  in  materials  science  (e.g.  the  Hall-Petch  effect  [74]).  In  the  analysis  that  follows 
next,  I  —  Inm  is  chosen,  consistent  with  prior  phase  field  studies  of  structural  transfor¬ 
mations  [5,8,9,39].  This  length  corresponds  physically  to  the  distance  from  an  externally 
unstressed  fracture  surface  or  twin  boundary  within  which  atoms  in  a  crystal  deviate  from 
their  ideal  static  positions.  Also,  I  =  0.01  is  regarded  as  fixed  for  consistency  with  results 
of  the  study  of  a  generic  solid  considered  most  often  earlier  in  this  paper.  These  choices 
lead  to  a  representative  domain  size  of  L  =  0.1  p,m  for  a  material  element  containing 
a  single  localized  fracture  or  fully  twinned  zone.  The  present  analysis  precludes  plastic 
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Table  2.  Material  properties  (/  =  Inm,  /  =  0.01). 


Material 

MO  (GPa) 

T  [J/m2] 

r  [J/m2] 

R 

TO 

KO 

Tip  [nm] 

Generic 

100 

1 

1 

1 

0.5 

50 

0.1 

Sapphire  (R  twin) 

167 

6 

0.125 

0.021 

0.202 

73.8 

0.071 

Sapphire  (B  twin) 

167 

40 

0.745 

0.019 

0.635 

95.1 

0.274 

Calcite  (e"*“  twin) 

36.7 

0.347 

0.183 

0.527 

0.694 

98.3 

0.130 

Magnesium  (tensile) 

19.0 

1.5 

0.117 

0.078 

0.130 

16.6 

0.049 

Table  3.  Model  predictions  and  comparisons  for  real  crystals  (GPa). 


Material 

Pc 

Pp 

Ah,  frac 

Frac.  data  [ref] 

Pl 

Ah,  twin 

PSF,  twin 

Twin  data  [ref] 

Generic 

4.6 

2.0 

15.9 

_ 

2.0 

8.0 

20 

_ 

Sapphire  (R) 

14.5 

6.3 

26.6 

7.9  [50] 

0.91 

5.4 

3.5 

1.0  [50] 

Sapphire  (B) 

37.5 

16.3 

26.6 

- 

2.23 

16.9 

5.4 

4.0  [50] 

Calcite 

1.6 

0.71 

5.8 

- 

0.52 

4.1 

2.8 

10“3  [42] 

Magnesium 

2.5 

1.1 

3.0 

- 

0.30 

0.39 

4.8 

0.08  [76] 

deformation  in  the  form  of  conventional  slip  of  dislocations  that  might  occur  (e.g.  in  ductile 
metal  crystals);  to  address  this,  the  kinematic  framework  of  Section  2.1  should  be  extended 
to  account  for  dislocation  glide  and  stored  defect  content  [75],  as  should  the  governing 
equations  of  the  phase  field  theory  [37,69].  Although  dislocation  slip  and  deformation 
twinning  are  often  competing  mechanisms,  omission  of  plastic  flow  is  deemed  reasonable  in 
the  present  application  since  twinning  is  much  easier  than  slip  in  calcite  [42],  rhombohedral 
twinning  is  preferred  over  pyramidal  slip  in  sapphire  [50]  and  tensile  twinning  is  preferred 
over  pyramidal  slip  in  magnesium  [41]. 

Predictions  of  criteria  (4.1)-(4.7)  are  compared  in  Table  3.  Also  listed  are  data  from 
other  sources.  For  sapphire,  shear  stress  data  for  fracture  resistance  on  R  planes  and  twinning 
resistance  on  B  and  R  systems  are  obtained  from  nonlinear  thermoelastic  analysis  of  shear 
strengths  observed  in  shock  physics  experiments  [50].  Data  on  twinning  resistance  for 
calcite  is  obtained  from  direct  shear  experiments  [42];  however,  this  stress  as  reported  may 
underestimate  the  true  twinning  resistance  since  it  is  an  average  value  for  the  entire  specimen 
that  omits  the  effect  of  local  stress  concentrations  near  the  twin  nucleus.  Data  on  twinning 
resistance  for  magnesium  is  obtained  from  atomic  simulations  [76] .  First  examining  fracture 
results.  Table  3  shows  that  Pq  >  Pp  for  all  considered  twin  systems  of  all  materials,  and 
that  Pp  provides  closest  agreement  with  the  test  datum  on  sapphire,  with  Pc  and  Pth,  fac 
greatly  overestimating  fracture  strength.  Next  examining  twinning  results,  it  is  observed 
that  Pp  <  Pth,  twin  and  Pp  <  Psp,  twin  for  all  materials.  For  R  twinning  in  sapphire  and 
for  twinning  in  calcite  and  magnesium,  Pp  is  closer  than  the  other  two  criteria  to  the  test 
data;  for  B  twinning  in  sapphire,  Pp  and  Psp  twin  are  of  comparable  accuracy  relative  to  the 
datum. 
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5.  Conclusions 

A  continuum  phase  field  theory  for  fracture  and  deformation  twinning  has  been  developed. 
Key  features  of  the  theory  include  assignment  of  distinct  order  parameters  for  fracture  and 
twinning  behaviours,  treatment  of  finite  deformation  and  nonlinear  elastic  response,  and 
degradation  of  elastic  strain  energy  and  twin  boundary  energy  with  cumulative  damage. 
The  present  paper  reports  the  first  known  model  accounting  for  all  such  features.  Out¬ 
comes  of  the  model  have  been  analysed  for  an  element  of  material  undergoing  simple 
shear  deformation.  Analytical  and/or  numerical  solutions  have  been  obtained  for  scenarios 
involving  isolated  physical  mechanisms  depicted  by  the  model,  including  homogeneous 
elasticity,  homogeneous  damage,  localized  fracture  and  localized  twinning.  In  addition  to 
providing  immediate  insight  into  pure  shear  behaviour,  these  solutions  can  be  used  later  for 
verification  of  more  advanced  numerical  methods  to  eventually  be  applied  towards  more 
complex  boundary  value  problems.  Stable  domains  of  homogeneous  solutions  have  been 
identified,  extending  previous  derivations  that  only  addressed  linear  elastic  behaviour  and 
tensile  deformations.  Coupled  fracture  and  twinning  mechanisms  have  also  been  considered. 
Criteria  for  shear  stresses  leading  to  fracture  and  twinning  have  been  derived  and  evaluated 
for  several  real  materials:  sapphire,  calcite  and  magnesium.  These  criteria  have  been  shown 
to  usually  demonstrate  closer  agreement  with  available  test  data  than  do  simple  criteria 
based  on  the  concept  of  theoretical  strength. 
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